Hydrodynamic modeling of a hypothetical river diversion near empire, louisiana

marsh porosity option is unable to properly deal with anisotropic micro- topography which may favor flow in one direction over another. 77 A) ADCIRC B) RMA2 Figure 5.29 Baseline Water Surface Elevations 78 Figure 5.30 WSE Difference for Baseline Runs, ADCIRC Minus RMA2 Figure 5.31 Velocity Magnitude Difference for Baseline Runs, ADCIRC Minus RMA2 79 CHAPTER 6. SUMMARY AND RECOMMENDATIONS 6.1 Results Summary Results of the sensitivity analyses show that generally both the RMA2 and ADCIRC models respond similarly to increases in the eddy viscosity parameter. For both models it is shown that there is a lower limit for the prescribed eddy viscosity below which the models become unstable. The WSE solution only shows significant sensitivity to the eddy viscosity within the diversion channel, where larger eddy viscosity values produce steeper surface slopes. The velocity solution is more sensitive to eddy viscosity. Larger eddy viscosity values produce flatter cross-channel velocity gradients for all of the channels which are resolved in the domain. The effect is most apparent in that peak velocities may be significantly reduced by large values of eddy viscosity. Comparison to an empirically determined cross-channel velocity profile for the diversion channel suggests that even the lowest values of eddy viscosity which produce stable results under- predict peak velocities in the channels in the domain. Considering this, the best choice for an eddy viscosity value is probably the lowest value that will produce stable results. Particle tracking results show, for a given model, that the overall flow patterns are not influenced significantly by the eddy viscosity, nor is the median residence time of the diverted water. The bottom friction sensitivity tests show significant sensitivity of the WSE solution when Manning’s n in increased over the entire depth range or increased only for shallow depths. Generally increased bottom friction causes steeper WSE gradients throughout the domain and therefore higher WSE closer to the diversion channel. This effect is more pronounced when the bottom drag is increased over the whole range of 80 depth. Higher WSE leads to wetting of some of the higher areas in the domain. The resulting changes in WSE gradient and wetted area cause a slightly larger percentage of the diverted water to follow steeper gradients to the north and west and ultimately leave the domain through Lake Grande Ecaille. This shift in flow patterns can be seen in the particle tracking results. The particle tracking results also show that the median retention time of the diverted water is significantly increased by the Overall n Increase scenario. Beyond these generalizations, the sensitivities of the models the bottom friction parameters are more model specific. This is largely due to the different wet/dry algorithms, particularly that ADCIRC completely blocks flow through dry elements which can cause significant changes in flow path as areas of the mesh become wet and open to flow. 6.2 Recommendations for Future Work While the results presented in this thesis may aid future modeling efforts in calibration by showing the solution sensitivities to some of the input parameters, observed data must be collected in order to calibrate and validate future models. Proper calibration and validation will require dynamic simulations including forcing from tides, other inflows, winds, rainfall, and possibly other processes. For the RMA2 model, simulations should be run using more advanced turbulence closure schemes and model sensitivities to their respective parameters should be tested. Perhaps for the ADCIRC model a more advanced turbulence closure scheme could be implemented. In addition, a future ADCIRC model should consider the channelizing effect that ADCIRC’s wet/dry scheme produces. Carefully smoothing model geometry could be done to help smooth 81 wet/dry transitions or perhaps a new algorithm similar to RMA2’s marsh porosity, could be included in the code and applied. A 100,000 cfs diversion at Empire, LA will likely affect the flow in the river and also circulation throughout Barataria Bay. Considering this, and the fact that multiple diversions may interact with each other, future work should include building a high- resolution model which encompasses a larger portion of the Mississippi River delta and surrounding estuaries. Additionally, the modeling of sediment transport and deposition in order to simulate the potential land building abilities of large scale river diversion will require long-term (tens of years or longer) simulations. These considerations will require more computational power and therefore will probably require the use of a parallel model which can take advantage of high performance computing. Parallel ADCIRC is already available and has been widely used for simulating ocean tides, hurricane storm surge, and even a river diversion into the Maurepas Swamp. A parallel version of RMA2 has also been created with limited modifications to the serial code (Rao 2005). However parallel RMA2 has not been widely used or tested on truly large-scale problems. As water is diverted from the river and flows through a diversion structure and into the diversion channel it can be expected that complex 3-D flows will develop and vertical acceleration and vertical variation in velocity will be important. This 3-D flow may have a significant effect on hydrodynamics in the diversion and also on the sediment which the diversion captures. Modeling of this will require the use of a fully three dimensional model which could be interfaced with a larger scale 2-D model. The modeling effort should be extended to include sediment transport modeling so that long term land building abilities of proposed river diversions can be analyzed. 82 This could be done using existing sediment transport codes such as SED2D, the sediment module within TABS-MD. Another option for future sediment transport modeling would be to enhance the particle tracking model and apply it as a model for sediment transport/deposition. It could also be parallelized making simulations with large numbers of particles and long durations practical. Particles could be used to represent parcels of sediment with specific properties (mass, density, fall velocity, size etc…). Also turbulent mixing could be considered with random particle displacements and a vertical velocity profile could be assumed as in the GISPART model. Based on these considerations, sediment parcels could be tracked and their deposition modeled to give an idea or where sediment from a diversion might deposit. Results from such a sediment deposition model could be synergistically combined with results from water quality models, ecological models, and physical models to help reveal the potential land building benefits of various river diversion scenarios which in the end could aid in restoration of our rapidly disappearing coast. 83 REFERENCES Barbé, Donald E., Kevin Fagot, John A. McCorquodale. (2000). Effects on Dredging Due to Diversions From the Lower Mississippi River. Journal of Waterway, Port, Coastal, and Ocean Engineering, 126(3), 121-129. Blanton, Brian. O. (1995). DROG3D: User’s Manual for 3-Dimensional Drogue Tracking on a Finite Element Grid with Linear Finite Elements. Ocean Processes Numerical Methods Laboratory Curriculum In Marine Science University of North Carolina at Chapel Hill. Britsch, Louis D., and Joseph B. Dunbar. (1993). Land Loss Rates: Louisiana Coastal Plain. Journal of Coastal Research 9(2), 324-338. Capps, Stephen A. (2003). Modeling a Mississippi River Diversion Into a Louisiana Wetland. MS. Thesis, Louisiana State University, Baton Rouge, LA. Connor, J.J., and C.A. Brebia. (1976). Finite Element Techniques for Fluid Flow. Butterworth & Co. Ltd. Woburn, Massachusetts, USA. Day, John W., Jr., Didier Pont, Philippe F. Hensel, and Carles Ibanez. (1995). Impacts of Sea-Level Rise on Deltas in the Gulf of Mexico and the Mediterranean: The Importance of Pulsing Events to Sustainability. Estuaries, 18(4), 636-647. DeLaune, R. D., A. Jugsujinda, G. W. Peterson, and W. H. Patrick Jr. (2003). Estuarine, Coastal and Shelf Science, 58, 653-662. Donnell, Barabara, Joseph V. Letter, Jr, William H. McAnally, and William A. Thomas. (2005). Users Guide to RMA2 WES Version 4.5. U.S. Army, Engineering Research and Development Center Waterways Experiments Station. Edwards, K. P., J. A. Hare, F. E. Werner, B. O. Blanton. (2006). Lagrangian circulation on the Southeast U.S. Continental Shelf: implications for larval dispersion and retention, Continental Shelf Research, 26(12-13), 1375-1394. King, Ian P., and Lisa C. Roig. (1991). Finite Element Modeling of Flow in Wetlands. Hydraulic Engineering. Proceedings of the National Conference on Hydraulic Engineering, 286-291. Kesel, Richard H. (1988). The Decline in the Suspended Load of the Lower Mississippi River and its Influence on Adjacent Wetlands, U.S.A. Environmental Geology and Water Science, 11(3), 271-281 Kesel, Richard H. (1989). The Role of the Mississippi River in Wetland Loss in Southeastern Louisiana, U.S.A. Environmental Geology and Water Science, 13(3), 183- 193 84 Kesel, Richard H., Elaine G. Yodis, and David J. McCraw. (1992). An Approximation of the Sediment Budget of the Lower Mississippi River Prior to Major Human Modification. Earth Surface Processes and Landforms, 17, 711-722. King, Ian P., and Lisa C. Roig. (1996). The Use of an Equivalent Porosity Method to Model Flow in Marshes. North American Water and Environment Congress, 3734-3739. Kundu, Pijush K., and Ira M. Cohen. (2004). Fluid Mechanics, third edition. Elsevier Academic Press, San Deigo, California, USA. Lane, Robert R., John W. Day Jr., and Jason N. Day. (2006). Wetland Surface Elevation, Vertical Accretion, and Subsidence at Three Louisiana Estuaries Receiving Diverted Mississippi River Water. Wetlands, 26(4), 1130-1142. Letter, Joseph V. and Gregory H. Nail. (1997). Wetland Engineering in Coastal Louisiana: Naomi Siphon. U.S. Army Corps of Engineers, Wetlands Research Program Technical Note, WG-RS-7.2. Letter, Joseph V. (1997). Wetland Engineering in Coastal Louisiana: Mississippi River Delta Splays. U.S. Army Corps of Engineers, Wetlands Research Program Technical Note, WG-RS-7.1. Luettich, Rick, and Joannes Westerink. (2000). ADCIRC, A (Parallel) ADvanced CIRculation Model for Oceanic, Coastal and Estuarine Waters. U.S. Army, Engineering Research and Development Center Waterways Experiments Station. Luettich, Rick, and Joannes Westerink. (2004). Formulation and Numerical Implementation of the 2D/3D ADCIRC Finite Element Model Version 44.XX. Mashriqui, Hassan S., and G. Paul Kemp. (2004) Hydrodynamic Models of Subprovince 2. Louisiana Coastal Area(LCA) report, vol.4, app. C, ch. 4, c54-c85. Morris, James T., P. V. Sundareshwar, Christopher T. Nietch, Björn Kjerfve, and D. R. Canhoon. (2002). Responses of Coastal Wetlands to Rising Sea Level. Ecology, 83(10), 2869-2877. Park, D., M. Inoue, W. J. Wiseman, Jr., D. Justic, and G. Stone. 2004. High-Resolution Integrated Hydrology-Hydrodynamic Model: Development and Application to Barataria Basin, Louisiana. U. S. Dept. of the Interior, Minerals Management Service, Gulf of Mexico OCS Region, New Orleans, LA. OCS Study MMS2004 - 063. 74 pp. Periañez, R. (2005). GISPART: a numerical model to simulate the dispersion of contaminants in the strait of Gibraltar. Environmental Modeling & Software, 20, 797- 802. 85 Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,. (1992). Numerical Recipies in c: The Art of Scientific Computing Second Edition, Cambridge University Press. Cambridge, UK. Rao, Prasada. (2005). A Parallel RMA2 model for simulating large-scale free surface flows. Environmental Modeling & software, 20(1), 47-53. Reed, Denise J. (1989). Patterns of Sediment Deposition is Subsiding Coastal Salt Marshes, Terrebonne Bay, Louisiana: the Role of Winter Storms. Estuaries, 12(4), 222- 227. Reed, Denise J. (1992). Effects of Weirs on Sediment Deposition in Louisiana Coastal Marshes. Environmental Management, 16(1), 55-65. Roig, Lisa C. (1994). Hydrodynamic Modeling of Flows in Tidal Wetlands. PhD Dissertation. University of Californina, Davis, CA. Roig, Lisa C. (1995). Mathematical Theory and Numerical Methods for the Modeling of Wetland Hydraulics. Water Resources Engineering, Proceedings of the First International Conference, San Antonio, TX, ASCE, 249-253. Turner, Eugene R., Joseph J. Baustian, Erick M. Swenson, and Jennifer S. Spicer. (2006). Wetland Sedimentation from Hurricanes Katrina and Rita. Science, 314(5798), 449-452. URS Corporation. (2006). Diversion Modeling Report, Mississippi River Re-introduction Into Maurepas Swamp, PO-29, DRAFT, URS Corporation December. 2006. Wilkerson, Gregory V., Jessican L. McGahan. (2005). Depth-Averaged Velocity Distribution is Straight Trapezoidal Channels. Journal of Hydraulic Engineering Technical Note, 131(6), 509-511. Willson, Clinton S., Nathan Dill, William Barlett, Samantha Danchuk, and Ryan Waldron. (2007). Physical and Numerical Modeling of River and Sediment Diversion in the Lower Mississippi River Delta. ASCE Coastal Sediments 2007, 1, 749-761. 86 APPENDIX. PARTICLE TRACKING CODE // maureparticle // by Nathan Dill // 11/13/06 ////////////////////////////////////////////////// ////////////////////////////////////////////////// function maureparticle // this program tracks particles through a 2d unstructured triangular mesh // using a known time varying or steady velocity field. the velocity within // an element is interpolated linearly bwtween the three nodes of the element. // particles are "moved" to new positions using either a direct Eulerian step, // or using the 2nd order Runge-Kutta method. // // the mesh may be provided as an ADCIRC type mesh named fort.14, as a RMA2 // type mesh named mesh.geo, or as generic files nodes.txt and element.txt. // in the fort.14 case, "weirs" if present are converted to elements; in the // mesh.geo case quadrilateral elements are split into triangles. // in all cases a files newfort.14(viewable in SMS), and nodes.txt // and elements.txt are made representing the weirless trianguar mesh // used for tracking. // // element to element and node to element neighbor tables are created and // used to relocate particles when they move from one element to the next. // the creation of these tables can be time consuming for large meshes, but // it is only necessary to do this once for each mesh. particles that leave // the mesh domain are brought back to the boundary perpendicular to their // position and in this manner can be tracked along the boundary or may move // back into the domain if the velocity field permits. // // One must be careful to specify a tracking timestep small enough that // particles cannot move out of an entire node-element group in one step. // a general rule could be that the timestep should be short enough prevent // a particle from moving more than half way across an element in one step. // i.e. timestep = (shortest element side)/2/(highest speed) // if the timestep is too large particles may get "lost" // ///////////////////////////////////////////////// // files necessary for maureparticle: // // fort.14, mesh.geo, or nodes.txt and elements.txt // (nodes.txt and elements.txt are the files actually read and used for tracking) // // fort.64 (ADCIRC output) // or velocity.dat (velocity ASCII file generated by SMS from RMA2's *.sol file) // // particles.inp (ASCII file with particles starting position) // format is: xstart ystart locat (locat is the starting element index, // " " " if 0 is used for locat Maureparticle // " " " will find each particles starting // " " " element, but this may take some time) // ////////////////////////////////////////////////// // files created by maureparticle: // 87 // newfort.14 (this file is for viewing mesh in SMS, it is an ADCIRC style mesh // but does not contain any boundary information) // // nodes.txt(if other mesh is given initially), format is: nid x y z // " " " " // elements.txt(if other mesh is given initially), format is: N1 N2 N3 // " " " // position.out (this is the particle position output file) // format is: for i=1:stp // TS(i) // for k=1:np // x(i,k) y(i,k) locat(i,k) // end // end // position.out may be very large, so use daytracks function(included below) // to make a more manageable output file "daytracks.tab" // /////////////////////////////////////////////////// // How to run maureparticle: // use "getf maureparticle.sci" to load the function // then direct SCILAB to the directory in which the the input files reside // run maureparticle by typing "maureparticle" at the command line // then input the answers to the questions asked. the first time a particular // mesh is used it is necessary to make node and element tables, if the tables // have already been made you can skip this part. it is necessary to know // how may words(strings) are on the first line of fort.14, or on the first // three lines of mesh.geo, also for mesh.geo it is necessary to know the // total number of element is the mesh so that the files are read in properly // the progress of maureparticle can be monitored by watching // the position.out grow. If desired, once maureparticle is finished // run the daytracks function to make a more manageable output file. // stacksize(100000000) format('v',18); global x y z nn global N1 N2 N3 ne global el2el node2el global px py locat np global TS tts tot dyn velfil global vx vy vxp vyp ttime global vx1 vx2 vy1 vy2 deltat global midnods disp('need to make node and element tables for this grid?') mktbls=input('(yes=1,no=0)'); if mktbls==1 disp('what kind of mesh file?'); grdfil=input('nodes.txt&elements.txt=0; fort.14=1; mesh.geo=2'); end tts=input('tracking timestep?(seconds)'); tot=input('total tracking time(days)?'); tot=tot*24*3600; dyn=input('steady-state(0) or dynamic(1)?') velfil=input('fort.64(0) or velocity.dat(1)?') 88 if dyn==0 & velfil==0 TS=input('velocity solution timestep for SS tracking(seconds)') end method=input('eulerian(0) or RK2(1)?') if mktbls==1 // make el2el and node2el tables if grdfil==1 // ADCIRC type mesh, may have to convert weirs to elements removeweirs elseif grdfil==2 // RMA2 type mesh, may have to split quadrilateral elements splitquads else // generic mesh, read in node.txt file and elements.txt file maketables disp('writing newfort.14') fdd=mopen('newfort.14','w') frstline='newfort file' mfprintf(fdd,'%s\n',frstline) mfprintf(fdd,'%i %i\n',ne,nn) nid=[1:nn]'; eid=[1:ne]' mfprintf(fdd,'%i %18.8f %18.8f %18.8f\n',nid,x,y,z) three=zeros(ne,1); three(:,:)=3; mfprintf(fdd,'%i %i %i %i %i\n',eid,three,N1,N2,N3) mclose(fdd) clear('three','nid','eid') end else readfiles1 end readfiles2 /////////////////////////////////////////////////////////////////////////////////// // find starting elements for particles if locat(1)==0 disp 'locating particles' for i=1:np for j=1:ne // clculate cross product to find what element particle is in locat(i)=j; ds1=[x(N1(j))-px(i) y(N1(j))-py(i)]; ds2=[x(N2(j))-px(i) y(N2(j))-py(i)]; ds3=[x(N3(j))-px(i) y(N3(j))-py(i)]; cros1=det([ds1; ds2]); cros2=det([ds2; ds3]); cros3=det([ds3; ds1]); if cros1>=0 & cros2>=0 & cros3>=0 break end end end end 89 /////////////////////////////////////////////////////// // begin writing output out=mopen('position.out','w'); time=0; mfprintf(out,'%i\n',time); mfprintf(out,'%18.8f %18.8f %i\n',px,py,locat); //////////////////////////////////////////////////// // read inital velocity field if dyn==0 & velfil==0 disp('reading fort.64, steady-state') read64std end if dyn==0 & velfil==1 disp('reading velocity.dat, steady-state') readdatstd end if dyn==1 & velfil==0 disp('reading fort.64, dynamic') ttime=0; read64dyn end if dyn==1 & velfil==1 disp('reading velocity.dat, dynamic') ttime=0 readdatdyn end ///////////////////////////////////////////////////////// // track particles steady state if dyn==0 disp('tracking particles, steady-state') stp=tot/tts for iii=1:stp // track loop [vxp,vyp]=velinterp(px,py) if method==1 px1=px; py1=py; end // find new positions Euler method px=px+tts*vxp; py=py+tts*vyp; time=time+tts; if method==1 // RK2 method px2=px; py2=py; pxmid=(px1+px2)/2; pymid=(py1+py2)/2; [vxp,vyp]=velinterp(pxmid,pymid) px=px1+tts*vxp; py=py1+tts*vyp; elcheck // update elemental location, and 90 // if left boundary, the position else // using Euler method elcheck end mfprintf(out,'%i\n',time); mfprintf(out,'%18.8f %18.8f %i\n',px,py,locat); end // track loop, for iii end // if dyn==0 //////////////////////////////////////////////////////// //////////////////////////////////////////////////////// // track particles dynamic if dyn==1 disp('tracking particles, dynamic') stp=tot/tts for iii=1:stp // dynamic track loop if ttime==1 //read the next velocity solution timestep if necessary vx1=vx2; vy1=vy2; if velfil==0 [n,secs2,j1]=mfscanf(1,fd64,'%f%f'); [n,NID,vx2,vy2]=mfscanf(nn,fd64,'%i%f%f'); clear NID else [n,junk]=mfscanf(2,veldat,'%s'); clear junk; [n,hours2]=mfscanf(1,veldat,'%f') [n,one]=mfscanf(lne,veldat,'%i'); clear one; [n,vx2,vy2]=mfscanf(lnn,veldat,'%f%f') vx2(midnods)=[] vy2(midnods)=[] end ttime=0 end vx=vx1+(vx2-vx1)*ttime; vy=vy1+(vy2-vy1)*ttime; [vxp,vyp]=velinterp(px,py) if method==1 px1=px; py1=py; end // find new positions Euler method px=px+tts*vxp; py=py+tts*vyp; time=time+tts; if method==1 // RK2 method px2=px; py2=py; pxmid=(px1+px2)/2; 91 pymid=(py1+py2)/2; [vxp,vyp]=velinterp(pxmid,pymid) px=px1+tts*vxp; py=py1+tts*vyp; elcheck else // using Euler method elcheck end mfprintf(out,'%i\n',time); mfprintf(out,'%18.8f %18.8f %i\n',px,py,locat); ttime=ttime+tts/deltat end // dyn track loop, for iii end // if dyn==1 ////////////////////////////////////////////////////////// mclose('all') endfunction // maureparticle ///////////////////////////////////////////////////// //------------------------------------------------------------------------------------- ////////////////////////////////////////////////////// ///////////////////////////////////////////////////// function removeweirs // // this function converts any weirs in the ADCIRC mesh // into elements so that particles can be tracked over // the weirs if necessary. it also builds the el2el // and node2el tables. it writes out files: // el2el.tbl, node2el.tbl, nodes.txt, elements.txt // and newfort.14(without any boundaries) // // stacksize(10000000) format('v',18); //read first two lines words=input('how many words are on first line of fort.14?') fd=mopen('fort.14','r'); fd2=mopen('newfort.14','w'); ffd=mopen('nodes.txt','w'); fffd=mopen('elements.txt','w') [n,s1]=mfscanf(words,fd,'%s'); // see how many words are on first line of fort.14 mfprintf(fd2,'%s ',s1); // you may need to read more strings mfprintf(fd2,'\n') [n,ne,nn]=mfscanf(1,fd,'%f %f'); // read nodes disp('reading nodes') x=zeros(nn,1); y=zeros(nn,1); z=zeros(nn,1); for i=1:nn 92 [n,id,xx,yy,zz]=mfscanf(1,fd,'%i%s%s%s'); // read as strings so not to loose precision! x(i)=evstr(xx); y(i)=evstr(yy); z(i)=evstr(zz); nid(i)=id; end mfprintf(ffd,'%i %18.8f %18.8f %18.8f\n',nid,x,y,z) mclose(ffd) //read elements disp('reading elements') [n,eid,three,N1,N2,N3]=mfscanf(ne,fd,'%i%i%i%i%i'); elements=[]; elements=[N1, N2, N3]; //read open boundries [n,nopen,j1,j2,j3,j4,j5]=mfscanf(1,fd,'%i%s%s%s%s%s'); //read strings for " = number of open boundaries" [n,nopenn,j1,j2,j3,j4,j5,j6,j7]=mfscanf(1,fd,'%i%s%s%s%s%s%s%s'); for i=1:nopen [n,nobn,j1,j2,j3,j4,j5,j6,j7,lb]=mfscanf(1,fd,'%i%s%s%s%s%s%s%s%i'); [n,nobid]=mfscanf(nobn,fd,'%i'); end // read land boundaries create new elements from weirs [n,nland,j1,j2,j3,j4,j5]=mfscanf(1,fd,'%i%s%s%s%s%s'); [n,nlandn,j1,j2,j3,j4,j5,j6,j7]=mfscanf(1,fd,'%i%s%s%s%s%s%s%s'); cnt=0 for q=1:nland [n,nnp,typ]=mfscanf(1,fd,'%i%i'); if typ==24 | typ==4 [n,s1,s2,s3,s4,s5,s6,s7,s8,s9,s10]=mfscanf(1,fd,'%s%s%s%s%s%s%s%s%s%s'); [n,N1,N2,elev,Wc1,Wc2]=mfscanf(nnp,fd,'%i%i%f%f%f'); //create elements from nodestring pairs cnt=cnt+1 //disp('creating new elements from weir nodestring pairs') z1=z(N1); z2=z(N2); x1=x(N1); x2=x(N2); y1=y(N1); y2=y(N2); ds1=[(x1(2)-x1(1)),(y1(2)-y1(1))]; //determine which way is counter clockwise ds2=[(x2(1)-x1(1)),(y2(1)-y1(1))]; cross=ds1(1)*ds2(2)-ds1(2)*ds2(1); if cross>0 // N2 is on N1's left EE=[]; ee=[]; for i=1:nnp-1 93 ee=[N1(i) N1(i+1) N2(i); N2(i) N1(i+1) N2(i+1)]; EE=[EE;ee]; end else // N2 on right EE=[]; ee=[]; for i=1:nnp-1 ee=[N2(i) N2(i+1) N1(i); N1(i) N2(i+1) N1(i+1)]; EE=[EE;ee]; end end elements=[elements;EE]; elseif typ==3 | typ==13 | typ==23 [n,j1,j2,j3,j4,j5,j6,j7,i1]=mfscanf(1,fd,'%s%s%s%s%s%s%s%s'); [n,NID1,elev,Wc1]=mfscanf(nnp,fd,'%i%f%f'); else [n,j1,j2,j3,j4,j5,j6,j7,i1]=mfscanf(1,fd,'%s%s%s%s%s%s%s%s'); [n,NID1]=mfscanf(nnp,fd,'%i'); end end if cnt>0 disp(cnt,'created elements from nodestring pairs') end mclose(fd); ne=length(elements(:,1)); mfprintf(fd2,'%i %i\n',[ne nn]); mfprintf(fd2,'%i %18.8f %18.8f %18.8f\n',nid,x,y,z); eid=[1:ne]'; three=zeros(ne,1); three(:,:)=3; mfprintf(fd2,'%i %i %i %i %i\n',[eid three elements]); mclose(fd2); mfprintf(fffd,'%i %i %i\n',[elements]) mclose(fffd) maketables // make el2el and node2el endfunction //removeweirs ///////////////////////////////////////////////////// //--------------------------------------------------------------------------------------- ///////////////////////////////////////////////////// ///////////////////////////////////////////////////// function splitquads // // Nathan Dill 94 // 9/11/06 // this script reads the RMA2 *.GEO file, splits any // quadrilateral elements into two triangle elements, // then writes a files new_nodes.txt and new_elements.txt // with the new elements added at the end // also writes out newfort.14 for viewing in sms // and corner_index.txt which is an index // to the original corner node numbers necessary for // getting the proper velocities at corner // nodes only from velocity.dat format('v',18); stacksize(100000000); fd=mopen('mesh.geo','r'); //words=6; //ne=6776 words=input('how many words are on first lines of mesh.geo?') ne=input('how many elements in mesh?') [n,j1]=mfscanf(words,fd,'%s'); clear('j1'); // read elements disp('reading elements') [n,GE,eid,n1,n2,n3,n4,n5,n6,n7,n8,one,zero]=mfscanf(ne,fd,'%s%i%i%i%i%i%i%i%i%i%i%f'); clear('GE','one','zero'); // read nodes disp('reading nodes') [n,GNN,nid,xx,yy,zz]=mfscanf(-1,fd,'%s%i%s%s%s'); nn=length(nid) mclose(fd) x=zeros(nn,1); y=zeros(nn,1); z=zeros(nn,1); for i=1:nn x(i)=evstr(xx(i)); y(i)=evstr(yy(i)); z(i)=evstr(zz(i)); end clear('xx','yy','zz','GNN') //corner nodes NN1=n1; NN2=n3; NN3=n5; N4=n7; // remove midside nodes midnodes=[n2; n4; n6; n8;]; count=0 for i=1:length(midnodes) k=find(midnodes(i)==nid); nid(k)=[]; x(k)=[]; y(k)=[]; z(k)=[]; 95 if length(k)>0 count=count+1; end end if count>0 disp(count,'removed midside nodes') end /// create index to renumber elements. index=zeros(max(nid),1); for i=1:length(index) k=find(nid==i); index(i)=k end // redefine elements with indexed node numbers N1=index(NN1); N2=index(NN2); N3=index(NN3); disp('writing nodes.txt') ffd=mopen('nodes.txt','w') mfprintf(ffd,'%i %18.8f %18.8f %18.8f\n',nid,x,y,z); mclose(ffd) // split quadrillateral elements // try not to make thin triangles disp('splitting quadrillateral elements') quadels=find(N4~=0); newels=[] for i=1:length(quadels) j=quadels(i) xx=[x(N1(j));x(N2(j));x(N3(j));x(index(N4(j)))]; yy=[y(N1(j));y(N2(j));y(N3(j));y(index(N4(j)))]; ds1=(xx(1)-xx(3))^2+(yy(1)-yy(3))^2; ds2=(xx(2)-xx(4))^2+(yy(2)-yy(4))^2; if ds1<ds2 // connect N1 and N3 to form new elements newels=[newels;N1(j),N3(j),index(N4(j));N1(j),N2(j),N3(j)]; else // connect N2 and N4 newels=[newels;N2(j),N3(j),index(N4(j));N1(j),N2(j),index(N4(j))]; end end // add newels to old triangle elements from GEO file // keeping only corner nodes N1(quadels)=[]; N2(quadels)=[]; N3(quadels)=[]; 96 oldels=[N1,N2,N3]; numnewels=(length(newels(:,1)))/2 newels=[oldels;newels]; // this is all the elements disp(numnewels,'split quadrillateral elements') disp('writing elements.txt') fd2=mopen('elements.txt','w'); mfprintf(fd2,'%i %i %i\n',newels) mclose(fd2) disp('writing newfort.14') fdd=mopen('newfort.14','w') nid=[1:1:length(x)]' z=(-1)*z ne=length(newels(:,1)) nn=length(x) mfprintf(fdd,'%s\n','no_boundary_grid') mfprintf(fdd,'%i %i\n',ne,nn); mfprintf(fdd,'%i %f %f %f\n',nid,x,y,z) three=zeros(ne,1) three(:,:)=3; eid=[1:1:ne]' mfprintf(fdd,'%i %i %i %i %i\n',[eid three newels]) mclose(fdd) fd=mopen('corner_index.txt','w') mfprintf(fd,'%i\n',index) mclose(fd) maketables // make el2el and node2el endfunction //splitquads ///////////////////////////////////////////////////////// //----------------------------------------------------------------------------------- /////////////////////////////////////////////////////// /////////////////////////////////////////////////////// function maketables // // Nathan Dill // 10/05/06 // this function reads the nodes.txt and elements.txt files // and makes the el2el.tbl and node2el.tbl. taken and modified // from original maureparticle_1 script // // global x y z nn global N1 N2 N3 ne global el2el node2el global px py locat np global TS tts tot dyn velfil global vx vy vxp vyp ttime global vx1 vx2 vy1 vy2 deltat readfiles1 97 elements=[N1 N2 N3]; /// make node to element table disp('making node to element table'); disp('this may take some time'); fd4=mopen('node2el.tbl','w'); node2el=zeros(nn,12); for i=1:nn k=find(elements(:,1)==i | elements(:,2)==i | elements(:,3)==i); node2el(i,1:length(k))=k; mfprintf(fd4,'%i %i %i %i %i %i %i %i %i %i %i %i\n',node2el(i,:)); end mclose(fd4); ////// use find command to make element to element table !this may take some time! disp('making element to element table') disp('this will take even longer') fd3=mopen('el2el.tbl','w'); el2el=zeros(ne,3); for i=1:ne kk=[]; for j=1:3 k=find(elements(i,j)==elements(:,1)); kk=[kk k]; k=find(elements(i,j)==elements(:,2)); kk=[kk k]; k=find(elements(i,j)==elements(:,3)); kk=[kk k]; end q=find(kk==i); kk(q)=[]; w=[]; for m=1:length(kk) n=find(kk==kk(m)); if length(n)==1 w=[w n]; end if length(n)==2 w=[w n(1)]; end end kk(w)=[]; el2el(i,1:length(kk))=kk; mfprintf(fd3,'%i %i %i\n',el2el(i,:)); end mclose(fd3) endfunction // maketables //////////////////////////////////////////////////////////// //--------------------------------------------------------- /////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////// function readfiles1 // 98 // this function reads in nodes.txt and element.txt // // // global x y z nn global N1 N2 N3 ne global el2el node2el global px py locat np global TS tts tot dyn velfil global vx vy vxp vyp ttime global vx1 vx2 vy1 vy2 deltat disp('reading nodes.txt') fd=mopen('nodes.txt','r'); [n,nid,xx,yy,zz]=mfscanf(-1,fd,'%i %s %s %s'); nn=length(nid) x=zeros(nn,1); y=zeros(nn,1); z=zeros(nn,1); for i=1:nn x(i)=evstr(xx(i)); y(i)=evstr(yy(i)); z(i)=evstr(zz(i)); end clear('xx','yy','zz'); mclose(fd); disp('reading elements.txt') fd=mopen('elements.txt','r'); [n,N1,N2,N3]=mfscanf(-1,fd,'%i %i %i'); mclose(fd); ne=length(N1); endfunction //readfiles1 /////////////////////////////////////////////////////////////// //------------------------------------------------------------ //////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////// function readfiles2 // // this function reads node2el.tbl, el2el.tbl, and particles.inp // global x y z nn global N1 N2 N3 ne global el2el node2el global px py locat np global TS tts tot dyn velfil global vx vy vxp vyp ttime global vx1 vx2 vy1 vy2 deltat // read el2el table 99 disp('reading el2el.tbl') fd1=mopen('el2el.tbl','r'); [n,E1,E2,E3]=mfscanf(ne,fd1,'%i%i%i'); el2el=[E1,E2,E3]; clear E1 E2 E3; mclose(fd1); // read node2el table disp('reading node2el.tbl') fd2=mopen('node2el.tbl','r'); [n,e1,e2,e3,e4,e5,e6,e7,e8,e9,e10,e11,e12]=mfscanf(nn,fd2,'%i%i%i%i%i%i%i%i%i%i%i%i'); node2el=[e1,e2,e3,e4,e5,e6,e7,e8,e9,e10,e11,e12]; clear e1 e2 e3 e4 e5 e6 e7 e8 e9 e10 e11 e12; mclose(fd2); // read particle input disp('reading starting particle positions'); fd3=mopen('particles.inp','r'); [n,pxx,pyy,locat]=mfscanf(-1,fd3,'%s%s%i'); px=evstr(pxx); py=evstr(pyy); np=length(locat); clear pxx pyy; mclose(fd3) endfunction readfiles2 ///////////////////////////////////////////////////////////// //---------------------------------------------------------- /////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// function read64std // global x y z nn global N1 N2 N3 ne global el2el node2el global px py locat np global TS tts tot dyn velfil global vx vy vxp vyp ttime global vx1 vx2 vy1 vy2 deltat fd4=mopen('fort.64','r'); [n,s1,s2,s3,s4,s5,s6]=mfscanf(1,fd4,'%s%s%s%s%s%s'); clear s1 s2 s3 s4 s5 s6 [n,nts,j2,j3,j4,j5]=mfscanf(1,fd4,'%f%f%f%f%f'); clear j2 j3 j4 j5 for i=1:nts [n,secs,j1]=mfscanf(1,fd4,'%f%f'); [n,NID,vx,vy]=mfscanf(nn,fd4,'%i%f%f'); clear NID if secs==TS break end end mclose(fd4); 100 endfunction //read64std ////////////////////////////////////////////////////////////// //----------------------------------------------------------- ///////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////// function readdatstd // global x y z nn global N1 N2 N3 ne global el2el node2el global px py locat np global TS tts tot dyn velfil global vx vy vxp vyp ttime global vx1 vx2 vy1 vy2 deltat fd=mopen('corner_index.txt','r') [n,index]=mfscanf(-1,fd,'%i') mclose(fd) fd=mopen('velocity.dat','r') [n,junk]=mfscanf(5,fd,'%s') [n,lnn,junk,lne]=mfscanf(1,fd,'%i%s%i') [n,junk]=mfscanf(5,fd,'%s') clear junk [n,one]=mfscanf(lne,fd,'%i') clear one [n,vx,vy]=mfscanf(lnn,fd,'%f%f') mclose(fd) k=find(index==0) vx(k)=[] vy(k)=[] endfunction // readdatstd /////////////////////////////////////////////////////////////// //------------------------------------------------------------ ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// function read64dyn // // this function reads the first two timesteps of the fort.64 // and begins the velocity as a linear funciton of time. global x y z nn global N1 N2 N3 ne global el2el node2el global px py locat np global TS tts tot dyn velfil global vx vy vxp vyp ttime global vx1 vx2 vy1 vy2 deltat fd64=mopen('fort.64','r'); [n,s1,s2,s3,s4,s5,s6]=mfscanf(1,fd4,'%s%s%s%s%s%s'); clear s1 s2 s3 s4 s5 s6 101 [n,nts,j2,j3,j4,j5]=mfscanf(1,fd4,'%f%f%f%f%f'); clear j2 j3 j4 j5 [n,secs1,j1]=mfscanf(1,fd4,'%f%f'); [n,NID,vx1,vy1]=mfscanf(nn,fd4,'%i%f%f'); [n,secs2,j1]=mfscanf(1,fd4,'%f%f'); [n,NID,vx2,vy2]=mfscanf(nn,fd4,'%i%f%f'); clear NID deltat=secs2-secs1 endfunction // read64dyn //////////////////////////////////////////////////////////// //---------------------------------------------------------- ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// function readdatdyn // global x y z nn global N1 N2 N3 ne global el2el node2el global px py locat np global TS tts tot dyn velfil global vx vy vxp vyp ttime global vx1 vx2 vy1 vy2 deltat global midnods fd=mopen('corner_index.txt','r') [n,index]=mfscanf(-1,fd,'%i') mclose(fd) veldat=mopen('velocity.dat','r') [n,junk]=mfscanf(5,veldat,'%s') [n,lnn,junk,lne]=mfscanf(1,veldat,'%i%s%i') [n,junk]=mfscanf(4,veldat,'%s') [n,hours1]=mfscanf(1,veldat,'%f') clear junk [n,one]=mfscanf(lne,veldat,'%i') clear one [n,vx1,vy1]=mfscanf(lnn,veldat,'%f%f') midnods=find(index==0) vx1(midnods)=[] vy1(midnods)=[] [n,junk]=mfscanf(2,veldat,'%s'); clear junk; [n,hours2]=mfscanf(1,veldat,'%f') [n,one]=mfscanf(lne,veldat,'%i'); clear one; [n,vx2,vy2]=mfscanf(lnn,veldat,'%f%f') vx2(midnods)=[] vy2(midnods)=[] deltat=(hours2*3600)-(hours1*3600); endfunction // readdatdyn 102 //////////////////////////////////////////////////////////// //---------------------------------------------------------- /////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////// function [vxpp,vypp]=velinterp(ppx,ppy) // // this function lineraly interpolates the velocity at the current // position (px,py) of the particles within the element locat. // the components of velocity are considered separately as a plane // which passes through three points defined by the xNi,yNi,vjNi // of the element. each velocity component at the point px,py // is found as the elevation on the plane at that point. global x y z nn global N1 N2 N3 ne global el2el node2el global px py locat np global TS tts tot dyn velfil global vx vy vxp vyp ttime global vx1 vx2 vy1 vy2 deltat for i=1:np nod1=N1(locat(i)); nod2=N2(locat(i)); nod3=N3(locat(i)); // x-velocity AA=y(nod1)*(vx(nod2)-vx(nod3))+y(nod2)*(vx(nod3)-vx(nod1))+y(nod3)*(vx(nod1)-vx(nod2)); BB=vx(nod1)*(x(nod2)-x(nod3))+vx(nod2)*(x(nod3)-x(nod1))+vx(nod3)*(x(nod1)-x(nod2)); CC=x(nod1)*(y(nod2)-y(nod3))+x(nod2)*(y(nod3)-y(nod1))+x(nod3)*(y(nod1)-y(nod2)); DD=-AA*x(nod1)-BB*y(nod1)-CC*vx(nod1); vxpp(i)=-1*(AA*ppx(i)+BB*ppy(i)+DD)/CC; // y-velocity AA=y(nod1)*(vy(nod2)-vy(nod3))+y(nod2)*(vy(nod3)-vy(nod1))+y(nod3)*(vy(nod1)-vy(nod2)); BB=vy(nod1)*(x(nod2)-x(nod3))+vy(nod2)*(x(nod3)-x(nod1))+vy(nod3)*(x(nod1)-x(nod2)); CC=x(nod1)*(y(nod2)-y(nod3))+x(nod2)*(y(nod3)-y(nod1))+x(nod3)*(y(nod1)-y(nod2)); DD=-AA*x(nod1)-BB*y(nod1)-CC*vy(nod1); vypp(i)=-1*(AA*ppx(i)+BB*ppy(i)+DD)/CC; end endfunction // velinterp ////////////////////////////////////////////////////////////////////////////// //--------------------------------------------------------------------------- ///////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////// function elcheck // // this function checks to see if the particles are still in the same elements // if not if first refers to the el2el table, then to the node2el table // and it updates the locat. if particle has left a boundary it places it back // on the boundary. 103 global x y z nn global N1 N2 N3 ne global el2el node2el global px py locat np global TS tts tot dyn velfil global vx vy vxp vyp ttime global vx1 vx2 vy1 vy2 deltat for i=1:np j=locat(i); ds1=[x(N1(j))-px(i) y(N1(j))-py(i)]; ds2=[x(N2(j))-px(i) y(N2(j))-py(i)]; ds3=[x(N3(j))-px(i) y(N3(j))-py(i)]; cros1=det([ds1; ds2]); cros2=det([ds2; ds3]); cros3=det([ds3; ds1]); if cros1<=0 | cros2<=0 | cros3<=0 for q=1:3 // refer to element to element table to look for particle jjj=el2el(j,q); if jjj==0 break end ds1=[x(N1(jjj))-px(i) y(N1(jjj))-py(i)]; ds2=[x(N2(jjj))-px(i) y(N2(jjj))-py(i)]; ds3=[x(N3(jjj))-px(i) y(N3(jjj))-py(i)]; cros1=det([ds1; ds2]); cros2=det([ds2; ds3]); cros3=det([ds3; ds1]); if cros1>=0 & cros2>=0 & cros3>=0 //disp 'found el2el' j=jjj; locat(i)=j; break end end if locat(i)~=jjj // find the closest node and refer to node to element table j=locat(i); // and look in those elements ds1=[x(N1(j))-px(i) y(N1(j))-py(i)]; ds2=[x(N2(j))-px(i) y(N2(j))-py(i)]; ds3=[x(N3(j))-px(i) y(N3(j))-py(i)]; mag1=sqrt(ds1(1)^2+ds1(2)^2); mag2=sqrt(ds2(1)^2+ds2(2)^2); mag3=sqrt(ds3(1)^2+ds3(2)^2); mag=[mag1 mag2 mag3]; corners=[N1(j) N2(j) N3(j)]; k=find(mag==min(mag)); closest=corners(k); 104 search=node2el(closest,:); k=find(search==0);, search(k)=[]; k=find(search==j);, search(k)=[]; for q=1:length(search) j=search(q); ds1=[x(N1(j))-px(i) y(N1(j))-py(i)]; ds2=[x(N2(j))-px(i) y(N2(j))-py(i)]; ds3=[x(N3(j))-px(i) y(N3(j))-py(i)]; cros1=det([ds1; ds2]); cros2=det([ds2; ds3]); cros3=det([ds3; ds1]); if cros1>=0 & cros2>=0 & cros3>=0 //disp 'found node2el' locat(i)=j; break end end end end if locat(i)~=j // the particle has left the domain, change its position //disp 'lost particle, move to boundary' // it may have skipped over a node-element group j=locat(i); // in which case you need to make the timestep smaller ds1=[x(N1(j))-px(i) y(N1(j))-py(i)]; // displacement vectors to nodes of element j ds2=[x(N2(j))-px(i) y(N2(j))-py(i)]; ds3=[x(N3(j))-px(i) y(N3(j))-py(i)]; mag1=sqrt(ds1(1)^2+ds1(2)^2); mag2=sqrt(ds2(1)^2+ds2(2)^2); // find the nearest node and all elements around it mag3=sqrt(ds3(1)^2+ds3(2)^2); mag=[mag1 mag2 mag3]; corners=[N1(j) N2(j) N3(j)]; k=find(mag==min(mag)); closest=corners(k); corners(k)=[] search=node2el(closest,:); k=find(search==0);, search(k)=[]; ell=[search', el2el(search,:)]; k=find(ell(:,4)==0); // this should be the two elements on the boundary bel=ell(k,1); if length(bel)==1 // closest is not on boundary ki=1 xn1=x(corners(1)); yn1=y(corners(1)); // positions of closest boundary nodes xn2=x(corners(2)); yn2=y(corners(2)); thetas=atan((yn2-yn1),(xn2-xn1)); slope=tan(thetas); if slope==0 105 slope=0.0000000001 ; end // oops divide by zero is no no newx=(slope*(py(i)-yn1+slope*xn1)+px(i))/(slope^2+1); // intersection perpendicular to boundary segemet newy=py(i)-(newx-px(i))/slope; else nods1=[N1(bel(1)) N2(bel(1)) N3(bel(1))]; // these are the nodes of the bels nods2=[N1(bel(2)) N2(bel(2)) N3(bel(2))]; els1=el2el(bel(1),:); // these are the neighbors of the bels els2=el2el(bel(2),:); nodsels1=[N1(els1(1)) N2(els1(1)) N3(els1(1)); N1(els1(2)) N2(els1(2)) N3(els1(2))] nodsels2=[N1(els2(1)) N2(els2(1)) N3(els2(1)); N1(els2(2)) N2(els2(2)) N3(els2(2))] share1=[] for pp=1:3 // loop to find shared node k=find(nodsels1==nods1(pp)); if length(k)==2 share1=nods1(pp) end end k=find(nods1==closest); nods1(k)=[]; k=find(nods1==share1); nods1(k)=[]; share2=[] for pp=1:3 k=find(nodsels2==nods2(pp)); if length(k)==2 share2=nods2(pp) end end k=find(nods2==closest); nods2(k)=[]; k=find(nods2==share2); nods2(k)=[]; nxnods=[nods1 nods2]; // these are the boundary nodes on either side of closest ds1=[x(nods1)-px(i) y(nods1)-py(i)]; // these are displacement vectors to next nodes on boundary ds2=[x(nods2)-px(i) y(nods2)-py(i)]; mag=[sqrt(ds1(1)^2+ds1(2)^2) sqrt(ds2(1)^2+ds2(2)^2)]; ki=find(mag==min(mag)) nxnod=nxnods(ki); // ki is the index of bel nxnods(ki)=[] notnod=nxnods // the other node if ki==2 notki=1 else notki=2 end xn1=x(nxnod); yn1=y(nxnod); // positions of closest boundary nodes 106 xn2=x(closest); yn2=y(closest); thetas=atan((yn2-yn1),(xn2-xn1)); slope=tan(thetas); if slope==0 slope=0.0000000001 ; end // oops divide by zero is no no newx=(slope*(py(i)-yn1+slope*xn1)+px(i))/(slope^2+1); // intersection perpendicular to boundary segemet newy=py(i)-(newx-px(i))/slope; end // if length(bel)==1 // check to see if new position is actually between boundary nodes dsb=[xn1-xn2 yn1-yn2] // displacement between boundary nodes dsnx=[newx-xn1 newy-yn1] // displacment to nxnod magb=sqrt(dsb(1)^2+dsb(2)^2); magnx=sqrt(dsnx(1)^2+dsnx(2)^2); if magnx>magb xn1=x(notnod); yn1=y(notnod); // positions of closest boundary nodes xn2=x(closest); yn2=y(closest); thetas=atan((yn2-yn1),(xn2-xn1)); slope=tan(thetas); if slope==0 slope=0.0000000001 ; end // oops divide by zero is no no newx=(slope*(py(i)-yn1+slope*xn1)+px(i))/(slope^2+1); // intersection perpendicular to boundary segemet newy=py(i)-(newx-px(i))/slope; px(i)=newx;, py(i)=newy; // assign new particle position locat(i)=bel(notki) // and element location else px(i)=newx;, py(i)=newy; // assign new particle position locat(i)=bel(ki) //disp(locat(i)) // and new element location end end //if lost end endfunction // elcheck ///////////////////////////////////////////////////////////////////////// //--------------------------------------------------------------------- ////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////// function daytracks(tint) 107 // this function reads the position.out file generated by maureparticle // and makes a pid,x,y,t tab delimited file "dayTracks.tab" // this file is a much more manageable size and allows for easy // examination of individual particle tracks. // "tint" is the desired(less frequent) output interval in seconds // particle tracks are listed one after another. // time is converted from seconds to days // you must know the total tracking simulation time, the tracking timestep // and the number of particles. format('v',18); //stacksize(100000000); fd1=mopen('position.out','r'); tts=input('tracking timestep (seconds)?'); totaltime=input('total tracking time (seconds)?'); np=input('number of particles?')' numints=totaltime/tint fd3=mopen('dayTracks.tab','w'); j=0 kk=0 xx=zeros(np,numints) yy=zeros(np,numints) tt=zeros(np,numints) for i=1:(totaltime/tts) // build arrays for resorting particle // postion and time data. [n,T1]=mfscanf(1,fd1,'%f'); time=zeros(np,1); time(:,:)=T1; [n,x,y,locat]=mfscanf(np,fd1,'%s%s%s'); clear locat if(T1==j) kk=kk+1 j=j+tint x=evstr(x) y=evstr(y) xx(:,kk)=x yy(:,kk)=y tt(:,kk)=time end end tt=tt./24./3600 mfprintf(fd3,'%s\t%s\t%s\t%s\n','partID','x','y','days') 108 for ii=1:np // write out particle position particle after particle ID=zeros(numints,1) ID(:,:)=ii x=xx(ii,:)' y=yy(ii,:)' t=tt(ii,:)' mfprintf(fd3,'%i\t%f\t%f\t%f\n',ID,x,y,t) end mclose('all') // endfunction // daytracks 109 VITA Nathan Lamont Dill was born February 20, 1979, in Phildelphia, Pennsylvania, to Pamela Druzilla Zimmer Dill and Richard Lamont Dill. He lived in Willingboro, New Jersey, until age three when he moved to Williamsport, Pennsylvania, with parents and two sisters. He attended public school in Williamsport. In eighth grade he competed as a mathcounts mathlete placing 111 in the state, and he also earned a perfect score on the National Science Olympiad’s physical science test. During the summer before ninth grade some home-made rocket fuel for his model rockets accidentally set fire to his bedroom but he managed to put the fire out before any significant damage was done to the house. In 1997 he won a silver medal at the PIAA District IV track meet clearing a height of 13 ft. in the pole vault. He graduated from the Williamsport Area High School in June 1997. In the fall he matriculated at Bowdoin College in Brunswick, Maine, where he pursued a major in physics. In 2000 he was awarded NCAA Division III All New England in the pole vault. On October 6, 2001 he was married in Weymouth, Massachusetts, to Janice Louise Donovan of Weymouth, daughter of Martha Louise Ferguson Donovan and Richard Jerome Donovan. In May 2002 he received his Bachelor of Arts Degree from Bowdoin College. After completing his undergraduate schooling he worked for two years as a physics teacher at the Northfield Mount Hermon School in Northfield, Massachusetts. During the summer of 2004 he moved to Baton Rouge, Louisiana, and began graduate studies in civil engineering which he thought based on prior experience might be safer than rocket science. He currently resides in Baton Rouge with his wife and three sons; Aidan Lamont, Josiah Donovan, and Paytah Jerome.