PROGRAM MainPulse c kzktexas1.f c This version 28 October 1998 c Obtained from http://people.bu.edu/robinc/kzk c c Robin Cleveland c Department of Aerospace and Mechanical Engineering c Boston University c Boston, MA 02215 c robinc@bu.edu c c Mark Hamilton c Department of Mechanical Engineering c The University of Texas at Austin c Austin, TX c hamilton@mail.utexas.edu c c This file: kzktexas1.f c The fortran source code for a time-domain computer code developed at c The University of Texas at Austin to model axisymmetric sound beams in c fluids. The code is based on an augmented KZK equation that accounts for c nonlinearity, diffraction, thermoviscous absorption, and absorption and c dispersion due to an arbitrary number of relaxation phenomena. The c principal references for this code are the following: c c Y.-S. Lee, "Numerical solution of the KZK equation for pulsed finite c amplitude sound beams in thermoviscous fluids," Ph.D. dissertation, The c University of Texas at Austin (1993). c c Y.-S. Lee and M. F. Hamilton, "Time-domain modeling of pulsed c finite-amplitude sound beams," J. Acoust. Soc. Am., Vol. 97, pp. 906-917 c (1995). c c R. O. Cleveland, M. F. Hamilton, and D. T. Blackstock, "Time-domain c modeling of finite-amplitude sound in relaxing fluids," J. Acoust. Soc. c Am., Vol. 99, 3312-3318 (1996). c c The dissertation by Lee (1993) contains the original listing of the code, c description of algorithm, and test cases. Lee's dissertation is summarized c in the article by Lee and Hamilton (1995). The article by Cleveland et al. c (1996) describes the inclusion of relaxation and other modifications of the c code. c c Modifications: c c Added in code to account for absorption and dispersion due c to multiple relaxation processes. c c Change nonlinear subroutine to reduce the stringent step size c restriction in Lee's code. In addition the nonlinear subroutine c will now force the code to take a smaller step-size if the c waveform looks like it will shock up at the next step. Note also c the Nonlinear routine no longer needs the workspace arrays P1, P2, c and t. It now borrows the a and b arrays from Tridiag c c The lines that made the code land directly onto Target points have c been changed so that that step-size doesn't become too small. Mike c Averkiou has found that this can yield extremely small stepsizes c which lead to numerical errors. c c The code can now output data even if the Target point is in the c IBFD region of calculation. c c The tau array is now calculated in the main loop c c The input procedure has been somewhat changed. c 1 The code now lets the user turn effects on or off with flags c at the beginning of the input. c 2 The user must now enter the number of points across the face of c the piston. c 3 The user now specifies both IBFD and CNFD step sizes c 4 The initial waveform can now be specified by the user providing c a file which contains the waveform across the face of the piston. c The file must be of the form c Any number of lines of general comments. The @ sign in the c first column indicates numerical data is to start c @ c rho0 tau0 dtau numpts c P(tau0,rho0) c P(tau0+dtau,rho0) c ... c P(tau0+dtau*(numpts-1),rho0) c rho1 tau1 dtau numpts c ... c rhoN tauN dtau numpts c ... c P(tauN+dtau*numpts-1)) c ... c c where all variables are in dimensionless form: rho is the radial c location of the pressure waveform, tau the time of the first sample c in the waveform, dtau the sampling rate, numpts the number of samples c in the waveform, and P the pressure values. Note that the code c automatically takes care of the disortion of the time samples due c to the co-ordinate transformation, that is, for a flat piston c all the tauj terms will be the same. c If only one waveform is supplied the code assumes a uniform piston c and applies the waveform over the rest of the face of the piston c The input routine will linearly interpolate from the tau grid of c the input waveform to the tau grid of the programme, i.e., dtau of c the input file does not have to equal dtau of the programme. c The input algorithm does not interpolate in the rho direction. c It simply puts an input rho waveform into the closest programme rho c location. However, if the input waveforms are sparse, that is, there c are programme rho points that are skipped the code will interpolate c on to the missing rho points. The best strategy is to choose the c number of programme piston points to be a multiple of the number of c piston points in the input file. c Note that rho0 should be zero, i.e., on axis, and rhoN=1.0, i.e., c the edge of the piston. If the final rho is not equal to 1 then the c code assumes that the rest of the piston has the same amplituide c as the last rho point. This means to get a uniform piston you simply c have to give the code the on axis waveform and it will force the rest c of the piston face to have this waveform too. c Be warned that if the code is to use 100 rho points one must provide c 101 initial pressure waveforms over the rho points 0 through 100. c c The Q array that Lee use to calculate the CNFDdiffract routine was c found to be superfluous and has been removed - almost halving the c memory requirements of the code. In addition the sumF array of the c IBFDdiffract routine and the sumC array of the CNFDdiffract are never c used at the same time so now they share the same array space, called c sum in the routine Exec. implicit none integer outfile,maxTauPts,maxRhoPts,arraySize,MaxPtVal integer datafile, maxRanges, maxOffAxis integer maxRelax real fudge c These parameters are used to dimension the arrays used c throughout the program. (Note that MaxPtVal MUST be set c to the maximum of either MaxTauPts or MaxRhoPts). c The user, of course, when running the program c may specify that less points are actually used. c The fudge factor is for the accuracy in checking for c stopping points in sigma steps. c parameter (outfile=10) parameter (datafile=11) parameter (maxTauPts=9000) parameter (maxRhoPts=1600) parameter (maxPtVal=max(maxTauPts,maxRhoPts)) parameter (arraySize=(maxTauPts+1)*(maxRhoPts+1)) parameter (maxRanges=50) parameter (maxOffAxis=20) parameter (maxRelax=10) parameter (fudge=1e-4) c VARIABLE DEFINITIONS c The following is a fairly complete list of variable definitions. c note that most of these vars are only really used in Exec. c c nlflag = flag which is 1 if nonlinear effects are included c diflag = flag which is 1 if diffraction effects are included c tvflag = flag which is 1 if thermoviscous effects are included c mrflag = flag which is 1 if spreading effects are included c c expo = exponent in envelope function c cycles = envelope duration in cycles c Kabsorption = coefficient of absorption c Knonlinearity = coefficient of nonlinearity c numRelax = number of relaxation processes c thetanu = array of the theta's for the relaxation c Cnu = the Cnu constants for the relaxation processes c NumTauPts = number of divisions along tau axis c NumRhoPts = # of divisions along rho axis c pistonPts = # of radial (rho) divisions within c a piston radius c TauMin,TauMax = min/max values for tau (used to define c the "tau window") c RhoMax = max radial units c SigMax = maximum range sigma for which calculations c are desired c P = array holding REAL pressure vs. tau and rho c at a specific sigma c tau = array giving tau values c dSigma = step size in sigma direction c dSigIB = INITIAL sigma step, implicit backward method c dSigCN = INITIAL sigma step, Crank-Nicolson method c dRho = step size in rho direction c dTau = step size in tau direction c NumRanges = number of ranges at which beam patterns are desired c watchPoints = array containing points at which to output data. c USAGE: (i,0) ==> a sigma "watchpoint" c (i,j), j.NE.0 ==> a rho "watchpoint" AT c the sigma value c in (i,0) c OffAxisCounts = array containing number of rho watchpoints c for each sigma watchpoint. c e.g. Trans...nts(i) = 2 ===> c that there are 2 rho watches c for the ith sigma watchpoint. c ThisRangeID = index used to keep track of which patterns have c been provided already c ThisRange = the desired pattern range at ThisRangeID. c outputform = indicates whether the user has selected c automatic sidelobe data or custom data. c outFormat = allocated space to hold exec's formatting c IBFDsteps = total number of IBFD iterations to perform, c c Any other variables are explained as needed. Integer variables c such as i,j,k,n1,n2,n3 are used as indices and as counters where c appropriate. Note that the variables are declared HERE (in c MainPulse) only to reserve memory. How the memory is actually c used is determined by the subroutines and how the variables are c redefined there. c c The variables declared here are mainly for workspace and to c pre-allocate storage space. c integer nlflag, diflag, tvflag, mrflag integer NumTauPts, NumRhoPts, NumRanges, numRelax integer pistonPts, IBFDsteps real dSigIB, dSigCN real tauMin, dTau real Kabsorption, Knonlinearity real Cnu(1:maxRelax), thetanu(1:maxRelax) real SigMax, RhoMax, dRho real cycles integer mexpo integer outputform integer OffAxisCounts(1:maxRanges) integer rhoIndex(0:maxOffAxis) real WatchPoints(1:maxRanges,0:maxOffAxis) integer i c WORKSPACES c for Exec: real P(arraySize) real tau(0:maxTauPts) c for TriDiag: real a(0:maxPtVal),b(0:maxPtVal),c(0:maxPtVal) real soln(0:maxPtVal),rhs(0:maxPtVal) real beta(0:maxPtVal),gamma(0:maxPtVal) c for diffraction routines real sum(0:maxRhoPts) c for Nonlinear c real t(0:numTauPts) = a() c real Pold(0:numTauPts) = b() c c Input is used to obtain the sizes needed for the several arrays c used in KZKPULSE. c (BEGIN) Call Input (maxTauPts,maxRhoPts,maxRanges,maxOffAxis, maxRelax, + nlflag, diflag, tvflag, mrflag, datafile, + cycles,mexpo,Kabsorption,Knonlinearity, + numRelax, Cnu, thetanu, + tauMin,dTau, numTauPts, + pistonPts, numRhoPts,dRho, + dSigIB,dSigCN, sigMax, IBFDsteps, + numRanges, WatchPoints, + OffAxisCounts, outputForm,fudge) OPEN(unit=outfile,file='pulse.out',status='unknown') c Initialize tau array do 500 i=0,numTauPts tau(i) = tauMin + i*dTau 500 continue if (cycles.NE.0.0) then Call Initialize (P, tau, numTauPts, numRhoPts, pistonPts, + dRho, cycles, mexpo) else Call ReadInitial (P, tau, numTauPts, numRhoPts, pistonPts, + datafile) endif Call Exec (maxPtVal,maxRanges,maxOffAxis, + nlflag, diflag, tvflag, mrflag, + Kabsorption,Knonlinearity, + numRelax, Cnu, thetanu, + dTau, numTauPts, + pistonPts, numRhoPts,dRho, + dSigIB,dSigCN, sigMax, IBFDsteps, + numRanges, WatchPoints, + OffAxisCounts, rhoIndex, outputForm, + P, tau,sum, + a,b,c,rhs,soln,beta,gamma, outfile,fudge) CLOSE(outfile) stop end c c c SUBROUTINE Input (maxTauPts,maxRhoPts,maxRanges,maxOffAxis, + maxRelax,nlflag, diflag, tvflag, mrflag, + datafile, cycles, mexpo, Kabsorption, + Knonlinearity, numRelax, Cnu, thetanu, + tauMin,dTau, numTauPts, + pistonPts, numRhoPts,dRho, + dSigIB,dSigCN, sigMax, IBFDsteps, + numRanges, WatchPoints, + OffAxisCounts, outputForm,fudge) C C This routine accepts input from the standard input device. c Input is divided into sections: c source and environmental params c output grid c computational grid implicit none c subroutinue's input arguments integer maxTauPts, maxRhoPts, maxRanges, maxOffAxis, maxRelax c subroutine's output arguments integer nlflag, diflag, tvflag, mrflag real cycles integer mexpo, datafile integer numTauPts, numRhoPts, pistonPts, numRelax real Kabsorption, Knonlinearity, tauMin, dTau, dRho real Cnu(1:maxRelax), thetanu(1:maxRelax) real dSigIB, dSigCN, sigMax integer IBFDsteps integer numRanges real Watchpoints(1:maxRanges,0:maxOffAxis) integer OffAxisCounts(1:maxRanges) integer outputForm real fudge c internal variables character*40 datafilename real tauMax, tauWindowC, Goldberg real pi,rTmp1, rTmp2, sigma, rho, RhoMax, minRhoMax integer i,j,k,iTmp1, iTmp2 integer LobeOutput, CustomOutput parameter (LobeOutput=1, CustomOutput=2) C VARIABLES c MOST variables are already defined in the Main routine. c here are the new ones: c iTmp1,rTmp1 = integer/real "dummy" variables c tauMax = max tau for tau window c tauWindowC = width of tau window on cycles c minRhoMax = min allowable RhoMax c c format statements... 1000 format(1X,A,G12.4) 1002 format(1X,A,G12.4,/) 1005 format(1X,A,I6,/) 1007 format(1x,a,i5,1x,a) 1009 format(1x,i5,1x,a) c BEGIN pi = acos(-1.0) write(*,*) 'KZKPULSE Program' write(*,*) '#################' write(*,*) ' ' write(*,*) 'Where appropriate, variable names appear in ' write(*,*) 'parenthesis with recommended values.' write(*,*) ' ' write(*,*) ' ' C C control parameters C ^^^^^^^^^^^^^^^^^^ write(*,*) 'CONTROL PARAMETERS' write(*,*) '===================' write(*,*) 'Nonlinear effects (1=yes):' read(*,*) nlflag write(*,*) 'Diffraction effects (1=yes):' read(*,*) diflag write(*,*) 'Thermoviscous effects (1=yes):' read(*,*) tvflag write(*,*) 'Relaxation effects (1=yes):' read(*,*) mrflag C C source and environmental parameters C ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ write(*,*) 'ENVELOPE PARAMETERS' write(*,*) '===================' write(*,*) 'ENTER the number of cycles (n) in the envelope:' read(*,*) cycles if (cycles.NE.0.0) then write(*,*) 'ENTER the Exponent (m) for the envelope:' write(*,*) '((m) controls the number of "rise-time cycles"(N)' write(*,*) 'through the equation N = n/m )' read(*,*) mexpo write(*,1000) 'ENVELOPE CYCLE = ',cycles write(*,1005) 'M EXPONENT = ',mexpo write(*,1002) 'RISE-TIME CYCLES = ',(cycles/mexpo) else write(*,*) 'ENTER the file name of the input data:' read(*,'(A40)') datafilename write(*,*) 'Fname ',datafilename open(unit=datafile,file=datafilename,status='old') mexpo=1.0 endif write(*,*) 'ENVIRONMENTAL PARAMETERS' write(*,*) '========================' write(*,*) 'ENTER Absorption and nonlinearity coefficients: ' read(*,*) Kabsorption,Knonlinearity write(*,1000) 'Absorption = ',Kabsorption write(*,1002) 'Nonlinearity = ',Knonlinearity Goldberg = Knonlinearity/Kabsorption write(*,*) 'ENTER the number of relaxation processes: ' write(*,1007) '[max=',maxRelax,']' read(*,*) numRelax write(*,1007) 'ENTER Cnu and thetanu for the ', + numRelax, 'processes' write(*,*) 'Cnu thetanu' do 2010 i=1,numRelax read(*,*) Cnu(i), thetanu(i) 2010 continue C C output grid parameters C ^^^^^^^^^^^^^^^^^^^^^^ write(*,*) 'OUTPUT PARAMETERS' write(*,*) '=================' 1020 write(*,*) ' ' write(*,*) 'CHOOSE desired OffAxis output type:' write(*,*) ' 1 - automatic selection of axial lobe plus ' write(*,*) ' the first two sidelobes' write(*,*) ' 2 - custom entry of OffAxis points' read(*,*) outputform if ((outputform.LT.1).OR.(outputform.GT.2)) goto 1020 write(*,*) 'ENTER the number of axial ranges desired:' write(*,1007) '[max=',maxRanges,']' read(*,*) numRanges if (numRanges.LT.1) numRanges=1 if (numRanges.GT.maxRanges) numRanges=maxRanges if (outputform.EQ.LobeOutput) then write(*,1007) 'ENTER the ',numRanges, + 'sigma range(s), in ASCENDING order.' sigMax = 0.0 do 1030 i=1,numRanges 1025 read(*,*) sigma if (sigma.LT.0.0) then write(*,*) 'POSITIVE values only. RE-ENTER.' goto 1025 endif if ((i.NE.1).AND.(sigma.LT.watchpoints(i-1,0))) then write(*,*) 'Enter values in ASCENDING order. RE-ENTER:' goto 1025 endif watchpoints(i,0) = sigma if (sigma.GT.sigMax) sigMax = sigma c store the RHO values of sidelobes at this sigma OffAxisCounts(i) = 2 watchpoints(i,1) = 5.1356*sigma/(1.0+sigma)/2.0 watchpoints(i,2) = 8.4172*sigma/(1.0+sigma)/2.0 1030 continue elseif (outputform.EQ.CustomOutput) then write(*,*) 'On the following lines, you will enter' write(*,*) 'each desired SIGMA RANGE, followed by the NUMBER' write(*,*) 'of Off-Axis points you desire for that range.' write(*,*) 'Following that, you will be asked to input the' write(*,*) 'NORMALIZED (XI=r/a) transverse radial positions' write(*,*) 'for your desired OffAxis points.' write(*,*) ' ' write(*,*) 'Be sure to enter the sigma ranges in an' write(*,*) 'ASCENDING order. NOTE that a maximum of ' write(*,1009) maxOffAxis,'Off-Axis points are accepted,' write(*,*) 'and that On-axis data is AUTOMATICALLY included.' write(*,*) '(so you need not enter a Off-Axis coord of 0.)' write(*,*) ' ' sigMax = 0.0 do 1050 i=1,numRanges 1035 write(*,1007) 'Pair#',i, + ': Enter sigma and # of Off Axis points:' read(*,*) sigma, OffAxisCounts(i) if (sigma.LT.0.0) then write(*,*) 'POSITIVE values only. RE-ENTER.' goto 1035 endif if ((i.NE.1).AND.(sigma.LT.watchpoints(i-1,0))) then write(*,*) 'Enter values in ASCENDING order. RE-ENTER:' goto 1035 endif if (OffAxisCounts(i).LT.0) OffAxisCounts(i) = 0 if (OffAxisCounts(i).GT.maxOffAxis) + OffAxisCounts(i) = maxOffAxis watchpoints(i,0) = sigma if (sigma.GT.sigMax) sigMax = sigma c obtain the OffAxis points, if any are desired, and convert to RHO. if (OffAxisCounts(i).GT.0) then write(*,1007) 'Enter the ',OffAxisCounts(i), + 'Off Axis position(s), in XI = r/a coordinates' read(*,*) (WatchPoints(i,j),j=1,OffAxisCounts(i)) c filter out zero and negative values j=1 do while (j.LE.OffAxisCounts(i)) if (WatchPoints(i,j).LE.0.0) then do 1037 k=(j+1),OffAxisCounts(i) watchPoints(i,k-1) = watchPoints(i,k) 1037 continue OffAxisCounts(i) = OffAxisCounts(i) - 1 else j = j + 1 endif enddo c convert to rho, if any OffAxis are left after the above filtering if (OffAxisCounts(i).GT.0) then do 1040 j=1,OffAxisCounts(i) rho = watchpoints(i,j)/(1.0 + sigma) watchpoints(i,j) = rho 1040 continue endif endif write(*,1009) OffAxisCounts(i),'radial positions accepted.' 1050 continue endif c now we must determine the maximum rho that is specified c from this data and ask the operator to enter at least c this much later on. minRhoMax = 0.0 do 1060 i=1,numRanges if (OffAxisCounts(i).GT.0) then do 1055 j=1,OffAxisCounts(i) if (watchpoints(i,j).GT.minRhoMax) + minRhoMax = watchPoints(i,j) 1055 continue endif 1060 continue c adjust by small amount to guard against computational c roundoff-error when making comparisons against these values. do 1065 i=1,NumRanges WatchPoints(i,0) = WatchPoints(i,0) - fudge 1065 continue C C numerical grid parameters C ^^^^^^^^^^^^^^^^^^^^^^^^^ c the default value for the piston pTs is pistonPts= 34 c The dSigs used to be hard-wired values. The user may now specify values, c these are used as default values if the user enters zeros dSigIB = 1.0e-3 dSigCN = 3.5e-3 c ask for other numerical grid params: write(*,*) write(*,*) write(*,*) 'COMPUTATIONAL GRID PARAMETERS' write(*,*) '=============================' c compute recommended tau values (in cycles) rTmp1 = cycles*pi*(log(1000.0))**(1.0/(2.0*mexpo)) TauMin = (-(rTmp1 + 10.0*pi)) / (2.0*pi) TauWindowC = (2.0*rTmp1 + 30.0*pi) / (2.0*pi) TauMax = TauMin + TauWindowC write(*,*) 'TAU WINDOW PARAMETERS...' write(*,*) write(*,*) 'Recommended min,max tau edges are in cycles:' write(*,*) 'TauMin <= ',TauMin,' cycles' write(*,*) 'TauMax >= ',TauMax,' cycles' write(*,*) 'ENTER Min. and Max Tau edge values in CYCLES' read(*,*) TauMin,TauMax tauWindowC = tauMax - tauMin write(*,1000) 'TAU MIN (cycles) = ', TauMin write(*,1000) 'TAU MAX (cycles) = ', TauMax write(*,1002) ' Cycles spanned = ', tauWindowC TauMax = TauMax * 2.0 * Pi TauMin = TauMin * 2.0 * Pi c compute recommended TauPtsPerCycle if (Goldberg.GE.10.0) then iTmp1 = 240 elseif (Goldberg.GE.1.0) then iTmp1 = 120 else iTmp1 = 60 endif 1070 continue write(*,*) 'Maximum total tau points = ',maxTauPts write(*,*) 'Recommended minimum: tau points/cycle = ',iTmp1 write(*,*) 'ENTER number of tau points to use PER CYCLE: ' read(*,*) iTmp2 numTauPts = nint (iTmp2 * tauWindowC) if (float(numTauPts).LT.(iTmp2*tauWindowC)) + numTauPts = numTauPts + 1 write(*,1005) 'TOTAL Tau Divisions = ',NumTauPts if (NumTauPts.GT.MaxTauPts) then write(*,*) '*** ERROR *** Maximum limit exceeded.' goto 1070 endif dTau = (tauMax-tauMin)/float(numTauPts) write(*,*) 'PISTON/"RADIUS" PARAMETERS' write(*,*) '==========================' 1080 continue 1085 write(*,*) 'The largest rho coordinate determined' write(*,*) ' from the data you entered was: ',minRhoMax write(*,*) 'To avoid unwanted reflections affecting' write(*,*) ' your data you should enter a RhoMax value' write(*,*) 'significantly larger than the value above,' write(*,*) 'unless you are examining such behavior.' write(*,*) ' ' write(*,*) 'A typical good value is 12.0.' write(*,*) ' ' write(*,*) 'ENTER maximum RHO coordinate for mesh: (RhoMax)' read(*,*) RhoMax if (RhoMax.LE.minRhoMax) then write(*,*) 'You MUST enter more than',minRhoMax goto 1085 endif write(*,*) RhoMax write(*,*) 'ENTER the number of points across the piston' write(*,*) 'The TOTAL number of radial computation points' write(*,1007) 'cannot exceed ',MaxRhoPts,'.' write(*,*) '(the default value is 100) :' read(*,*) itmp1 if (itmp1.NE.0) pistonPts=itmp1 write(*,1005) 'Number of points per piston radii = ',pistonPts dRho = 1.0/float(pistonPts) numRhoPts = nint( RhoMax* float(pistonPts)) if (float(numRhoPts).LT.(RhoMax*float(pistonPts))) + numRhoPts = numRhoPts + 1 write(*,1007) 'At ',pistonPts,'per rho, ' write(*,1005) 'TOTAL Rho Divisions = ',numRhoPts if (numRhoPts.GT.maxRhoPts) then write(*,*) '*** ERROR *** Maximum limit exceeded.' goto 1080 endif if (Goldberg.LT.10.0) then itmp1 = 50 else itmp1 = 100 endif write(*,*) 'SIGMA PARAMETERS' write(*,*) '================' write(*,*) ' ' write(*,*) 'ENTER the step size for the IB method (default=', + dSigIB,')?' read(*,*) rtmp2 c A fix so that old input files which don't have dSigIB and dSigCN c will still make the code run. if (rtmp2.GT.1.0) then IBFDsteps=rtmp2 write(*,1005) '# Implicit steps = ',IBFDsteps return endif if (rtmp2.GT.0.0) dSigIB=rtmp2 write(*,*) ' ' write(*,*) 'ENTER the number of implicit steps to perform' write(*,*) ' before switching to the Crank-Nicolson method' write(*,1007) ' (Recommended minimum = ',itmp1,')' read(*,*) IBFDsteps write(*,1005) '# Implicit steps = ',IBFDsteps write(*,*) ' ' write(*,*) 'ENTER the step size for the CN method (default=', + dSigCN,')?' read(*,*) rtmp2 if (rtmp2.GT.0.0) dSigCN=rtmp2 return end c c c SUBROUTINE Initialize (P, tau, NumTauPts, NumRhoPts, pistonPts, + dRho, cycles, mexpo) c This subroutine initializes the P implicit none integer numTauPts, numRhoPts real P(0:numTauPts,0:numRhoPts), tau(0:numTauPts) real dRho, cycles integer mexpo integer pistonPts integer i,j real pi c VARIABLE DEFINITIONS: (see Main) pi = acos(-1.0) c initialize the P grid that covers the piston area. Note that c the initialization function is done in terms of the retarded time c tau at t=0. do 13510 i=0,numTauPts do 13500 j=0,pistonPts P(i,j) = exp(-((tau(i)+(j*dRho)**2)/(cycles*pi))**(2*mexpo)) + * sin(tau(i)+(j*dRho)**2) 13500 continue 13510 continue c initialize the REST of the pressure grid to zero do 13530 j=(pistonPts+1),numRhoPts do 13520 i=0,numTauPts P(i,j) = 0 13520 continue 13530 continue return end c c c SUBROUTINE ReadInitial(P, tau, numTauPts,numRhoPts,pistonPts, + datafile) c This subroutine initialises the P array using initial data provided c by the user. implicit none integer numTauPts, numRhoPts real P(0:numTauPts,0:numRhoPts), tau(0:numTauPts) integer datafile, pistonPts character dummyc integer i,j,k integer jrho, jrholast ,jdrho, numpts real rho, dp real tauin, dtauin, delta, Pnext, Plast, taunext, taulast c VARIABLE DEFINITIONS: (see Main) c Read the P grid that covers the piston area from the input file provided c by the user, and set the rest of the intial conditions to zero. c Clear out comments at front of input file read(datafile,'(A)') dummyc do 13600 while (dummyc.NE.'@') read(datafile,'(A)') dummyc 13600 continue c Make sure first waveform is on axis read(datafile,*) rho, tauin, dtauin, numpts tauin=tauin-rho*rho jrho=nint(rho*float(pistonPts)) if (jrho.NE.0) then write(*,*) "Error in input data file. First pressure signature" write(*,*) "must be on axis. Programme is going to assume that" write(*,*) "that the first signature is on axis" jrho=0 endif jrholast=-1 c DO until file is run out of input points do 13680 while (jrho.GE.0) c write(*,*) 'dtauin',dtauin i=0 c Fill array with zeros until input waveform starts do 13610 while (tau(i).LT.tauin) P(i,jrho)=0.0 i=i+1 13610 continue c write(*,*) 'i tau(i) tauin',i,tau(i),tauin c Read in waveform and linearly interpolate onto tau array read(datafile,*) Pnext taunext=tauin k=1 do 13640 while (k.LT.numpts) taulast=taunext taunext=tauin + k*dtauin Plast=Pnext read(datafile,*) Pnext k=k+1 c write(*,*) 'k taulast taunext',k,taulast,taunext do 13620 while (tau(i).GT.taunext) if (k.GT.numpts) goto 13645 taulast=taunext taunext=tauin + k*dtauin Plast=Pnext read(datafile,*) Pnext k=k+1 c write(*,*) 'k taulast taunext',k,taulast,taunext 13620 continue dp=(Pnext-Plast)/dtauin do 13630 while (tau(i).LE.taunext) P(i,jrho)=Plast+(tau(i)-taulast)*dp i=i+1 13630 continue c write(*,*) 'i tau(i) taulast',i,tau(i),taulast 13640 continue c Fill the final part of waveform with zeros 13645 do 13650 i=i,numtaupts P(i,jrho)=0.0 13650 continue jdrho=jrho-jrholast if (jdrho.GT.1) then c Interpolate between rho points if input file doesn't have c a small enough drho do 13670 i=0,numtaupts dp=(P(i,jrho)-P(i,jrholast))/float(jdrho) do 13660 j=1,jdrho-1 P(i,jrholast+j)=P(i,jrholast)+dp*float(j) 13660 continue 13670 continue elseif (jdrho.LT.0) then write(*,*) "ReadInitial ERROR" write(*,*) "Initial data file has signatures that are not" write(*,*) "going monotonically away from the axis." stop endif jrholast=jrho read(datafile,*,END=13690) rho, tauin, dtauin, numpts jrho=nint(rho*float(pistonPts)) goto 13675 c Rewind datafile if no points left to fill up transducer 13690 if (jrholast.GE.pistonPts) goto 13700 jrho=jrho+1 write(*,*) "ReadInitial: axial waveform copied for rhopt:",jrho rewind(datafile) read(datafile,'(A)') dummyc do 13695 while (dummyc.NE.'@') read(datafile,'(A)') dummyc 13695 continue read(datafile,*) rho, tauin, dtauin, numpts rho=jrho/float(pistonPts) 13675 tauin=tauin-rho*rho 13680 continue 13700 close(datafile) c initialize the REST of the pressure grid to zero do 13830 j=(jrholast+1),numRhoPts do 13820 i=0,numTauPts P(i,j) = 0 13820 continue 13830 continue return end c c c SUBROUTINE Exec (maxPtVal,maxRanges,maxOffAxis, + nlflag, diflag, tvflag, mrflag, + Kabsorption,Knonlinearity, + numRelax, Cnu, thetanu, + dTau, numTauPts, + pistonPts, numRhoPts,dRho, + dSigIB,dSigCN, sigMax, IBFDsteps, + numRanges, WatchPoints, + OffAxisCounts, rhoIndex, outputForm, + P, tau,sum, + a,b,c,rhs,soln,beta,gamma, outfile,fudge) c c This routine is responsible for the main execution of the program. c It controls the implicit backward and Crank-Nicolson c finite differencing methods. c implicit none c passed arguments integer maxPtVal, maxRanges, maxOffAxis integer nlflag, diflag, tvflag, mrflag real Kabsorption, Knonlinearity integer numRelax real Cnu(1:numRelax), thetanu(1:numRelax) real dTau, dRho, dSigIB, dSigCN, sigmax, dSigShock integer numTauPts, numRhoPts, pistonPts, IBFDsteps real watchpoints(1:maxRanges,0:maxOffAxis) integer OffAxisCounts(1:maxRanges), outputform integer rhoIndex(0:maxOffAxis) integer numRanges integer LobeOutput, CustomOutput, outfile parameter (LobeOutput=1, CustomOutput=2) real fudge real P(0:NumTauPts,0:NumRhoPts) real tau(0:NumTauPts) c internal vars. real sigma real dSigma real Target integer TargetID integer i,j,iTmp, n1,n2 c sub-WORKSPACES: c for TriDiag: real a(0:maxPtVal),b(0:maxPtVal),c(0:maxPtVal) real soln(0:maxPtVal),rhs(0:maxPtVal) real beta(0:maxPtVal),gamma(0:maxPtVal) c for diffraction routines real sum(0:numRhoPts) c for Nonlinear c real t(0:numTauPts) = a() c real Pold(0:numTauPts) = b() c variable definitions for internal vars. c sigma = current sigma value c dSigma = sigma step size c Target = next sigma targetted for output by user c TargetID = index of the next target c in the watchpoints array c c (BEGIN) c prep the outfile with some data: write(outfile,*) 'PULSE OUTPUT FILE' write(outfile,*) 'Dimensional Data follows:' write(outfile,*) numRanges,numTauPts write(outfile,*) ((WatchPoints(n1,0)+fudge),n1=1,numRanges) write(outfile,*) 'RESULTS BEGIN:' write(outfile,*) ' ' c begin the math... write(*,*) ' ' write(*,*) 'Max Sigma = ',SigMax write(*,*) 'Beginning IBFD iterations...' Sigma = 0 TargetID=1 Target=WatchPoints(TargetID,0) c ================================================================= c implicit backward finite difference method is used c for diffraction and absorption. The nonlinearity is included c analytically. We perform IBFDsteps. c ================================================================= dSigShock=11*dSigIB do 10200 n1=1,IBFDsteps c Make sure step size is less than half plane wave shock formation c distance dSigma = min(((1.0+sigma)**2)*dSigIB,dSigShock/2.0) c Adjust step to fall on Target - but don't let dSigma get too small if ((Target-sigma).LT.dSigma) then dSigma = Target - sigma + fudge if (dSigma.LT.fudge) dSigma=fudge endif sigma = sigma + dSigma if (diflag.EQ.1) + call IBFDdiffract (P, numTauPts, numRhoPts, sigma, + dSigma, dTau, dRho, sum, + a,b,c,rhs,soln,beta,gamma) if (tvflag.EQ.1) + call IBFDabsorb (P, numTauPts, numRhoPts, dSigma, dTau, + Kabsorption, a,b,c,rhs,soln,beta,gamma) if (mrflag.EQ.1) then do 10150 n2=1,numRelax call IBFDrelax (P, numTauPts, numRhoPts, dSigma, dTau, + Cnu(n2), thetanu(n2), + a,b,c,rhs,soln,beta,gamma) 10150 continue endif if (nlflag.EQ.1) + call Nonlinear (P, sigma, dSigma, dTau, tau, a, b, + Knonlinearity, numTauPts, numRhoPts, + outfile, dSigShock) C SEND OUT THE OUTPUT DATA, ADJUST Target IF NECESSARY if (sigma.GE.Target) then write(outfile,*) 'PRESSURE DATA FOR SIGMA = ' write(outfile,*) sigma write(outfile,*) 'NUMBER OF RADIAL POINTS PROVIDED = ' write(outfile,*) OffAxisCounts(TargetID)+1 if (outputform.EQ.LobeOutput) then write(outfile,*) 'PRESSURE DATA ARE GIVEN FOR THE PEAKS' write(outfile,*) ' OF THE AXIAL & 1st & 2nd SIDELOBES.' else write(outfile,*) 'PRESSURE DATA ARE GIVEN ON-AXIS AND ' write(outfile,*) ' AT THE REQUESTED RADIAL POSITIONS' endif write(outfile,11050) rhoindex(0) = 0 if (OffAxisCounts(targetID).GT.0) then do 10180 j=1,OffAxisCounts(targetID) rhoIndex(j) = nint(watchpoints(targetID,j)*pistonPts) 10180 continue endif write(outfile,*) 0.0,(float(rhoindex(j))/float(pistonPts), + j=0,OffAxisCounts(targetID)) write(outfile,*) 0.0,((float(rhoindex(j))/float(pistonPts))* + (1.0+sigma), j=0,OffAxisCounts(targetID)) write(outfile,12010) do 10190 i=0,numTauPts write(outfile,*) tau(i),((P(i,rhoindex(j))/(1.0+sigma)), + j=0,OffAxisCounts(targetID)) 10190 continue write(outfile,*) ' ' write(*,*) 'Data stored for Sigma = ',sigma,' dS=',dSigma TargetID = TargetID + 1 Target = watchpoints(TargetID,0) if (sigma.GE.sigMax) go to 12000 endif 10200 continue write(*,*) 'Ended IBFD at sigma = ',sigma c ================================================================= c CRANK-NICOLSON finite difference method is used for diffraction c and absorption. The nonlinearity is included analytically c Output is only written for the requested ranges. c ================================================================= c If dSigma is too big, such that we'll leap over our next c desired range point, we'll adjust dSigma so that we c instead land on it, remembering to reverse the fudge bias c we introduced in the input routine. write(*,*) 'Beginning CNFD iterations...' C {REPEAT} c Make sure dsigma is less that half a shock formation distance 11000 dSigma = min(((1.0+sigma)**2)*dSigCN,dSigShock/2.0) cAdjust step to fall on Target - but don't let dSigma get too small if ((Target-sigma).LT.dSigma) then dSigma = Target - sigma + fudge if (dSigma.LT.fudge) dSigma=fudge endif sigma = sigma + dSigma if (diflag.EQ.1) + call CNFDdiffract(P, numTauPts, numRhoPts, sigma, + dSigma, dTau, dRho, sum, + a,b,c,rhs,soln,beta,gamma) if (tvflag.EQ.1) + call CNFDabsorb (P, numTauPts, numRhoPts, dSigma, dTau, + Kabsorption, a,b,c,rhs,soln,beta,gamma) if (mrflag.EQ.1) then do 11010 n2=1,numRelax call CNFDrelax (P, numTauPts, numRhoPts, dSigma, dTau, + Cnu(n2), thetanu(n2), + a,b,c,rhs,soln,beta,gamma) 11010 continue endif if (nlflag.EQ.1) + call Nonlinear (P, sigma, dSigma, dTau, tau, a, b, + Knonlinearity, numTauPts, numRhoPts, + outfile, dSigShock) C SEND OUT THE OUTPUT DATA, ADJUST Target IF NECESSARY if (sigma.GE.Target) then write(outfile,*) 'PRESSURE DATA FOR SIGMA = ' write(outfile,*) sigma write(outfile,*) 'NUMBER OF RADIAL POINTS PROVIDED = ' write(outfile,*) OffAxisCounts(TargetID)+1 if (outputform.EQ.LobeOutput) then write(outfile,*) 'PRESSURE DATA ARE GIVEN FOR THE PEAKS' write(outfile,*) ' OF THE AXIAL & 1st & 2nd SIDELOBES.' else write(outfile,*) 'PRESSURE DATA ARE GIVEN ON-AXIS AND ' write(outfile,*) ' AT THE REQUESTED RADIAL POSITIONS' endif write(outfile,11050) 11050 format(1x,/,1x,'TAU',T20,'Actual RHO',/,T20,'and XI [r/a]') rhoindex(0) = 0 if (OffAxisCounts(targetID).GT.0) then do 11100 j=1,OffAxisCounts(targetID) rhoIndex(j) = nint(watchpoints(targetID,j)*pistonPts) 11100 continue endif write(outfile,*) 0.0,(float(rhoindex(j))/float(pistonPts), + j=0,OffAxisCounts(targetID)) write(outfile,*) 0.0,((float(rhoindex(j))/float(pistonPts))* + (1.0+sigma), j=0,OffAxisCounts(targetID)) write(outfile,12010) 12010 format(1x,75('=')) do 11200 i=0,numTauPts write(outfile,*) tau(i),((P(i,rhoindex(j))/(1.0+sigma)), + j=0,OffAxisCounts(targetID)) 11200 continue write(outfile,*) ' ' write(*,*) 'Data stored for Sigma = ',sigma,' dS=',dSigma TargetID = TargetID + 1 Target = watchpoints(TargetID,0) endif C {UNTIL SIGMA >= SIGMAX} if (sigma.LT.sigMax) go to 11000 12000 return end c c c SUBROUTINE IBFDabsorb (P, numTauPts, numRhoPts, dSigma, dTau, + Kabsorption, a,b,c,rhs,soln,beta,gamma) c c This routine performs the implicit backward finite difference c method for the absorption. c c The routine assumes that each element in the array affects c only those elements adjacent to it in retarded time tau, NOT c those adjacent to it radially (in rho). Therefore, each c "rho line" in the array will only affect the corresponding c rho line at the next sigma step, not those adjacent to it. c c The rho lines thus become independent of each other, and allow the c finite difference to be carried out on each row independently. c Therefore, to conserve memory usage, after computing the solutions c for a particular rho line, the solutions are stored back into the c original P array, taking advantage of the rho line's independence. c implicit none integer numTauPts, numRhoPts real P(0:numTauPts,0:numRhoPts) real dSigma, dTau, Kabsorption real a(0:numTauPts),b(0:numTauPts),c(0:numTauPts) real rhs(0:numTauPts), soln(0:numTauPts) real beta(0:numTauPts), gamma(0:numTauPts) integer start, finish real lambda integer i,j c VARIABLE DEFINITIONS: See Main c and c lambda = a convenient constant relating dSigma and dTau c a,b,c,rhs,soln = temporary workspace used for call to Tridiag c C (BEGIN) start = 1 finish = numTauPts - 1 c initialize the coefficient "bands" in the tridiagonal matrix. c these "bands" remain the same regardless. lambda = Kabsorption*dSigma/(dTau**2) do 20010 i=start,finish a(i) = -lambda b(i) = 1.0 + 2.0*lambda c(i) = -lambda 20010 continue c Here we will solve the tridiagonal matrix for each "rho" coordinate, c thereby generating solutions for the P values with the same "rho" c but at the next sigma step. First the right-hand-side matrix is c initialized, then the tridiagonal routine is called to solve it, and c solutions are stored back into the P array to record the P values c at the NEXT sigma step. Each "rho" line is assumed to be independent c of the other rho lines, which allows us to solve this as a series of c 2-dimensional finite difference systems (sigma, tau). cvd$ cncall do 20040 j=0,(numRhoPts-1) do 20020 i=start,finish rhs(i) = P(i,j) 20020 continue call Tridiag(start,finish,a,b,c,rhs,soln,beta,gamma,numTauPts) do 20030 i=start,finish P(i,j) = soln(i) 20030 continue 20040 continue return end c c c SUBROUTINE IBFDrelax (P, numTauPts, numRhoPts, dSigma, dTau, + Cnu,thetanu ,a,b,c,rhs,soln,beta,gamma) c c This routine performs the implicit backward finite difference c method for a single relaxation process. c c The routine assumes that each element in the array affects c only those elements adjacent to it in retarded time tau, NOT c those adjacent to it radially (in rho). Therefore, each c "rho line" in the array will only affect the corresponding c rho line at the next sigma step, not those adjacent to it. c c The rho lines thus become independent of each other, and allow the c finite difference to be carried out on each row independently. c Therefore, to conserve memory usage, after computing the solutions c for a particular rho line, the solutions are stored back into the c original P array, taking advantage of the rho line's independence. c implicit none integer numTauPts, numRhoPts real P(0:numTauPts,0:numRhoPts) real dSigma, dTau, Cnu, thetanu real a(0:numTauPts),b(0:numTauPts),c(0:numTauPts) real rhs(0:numTauPts), soln(0:numTauPts) real beta(0:numTauPts), gamma(0:numTauPts) integer start, finish real lambda, mu integer i,j c VARIABLE DEFINITIONS: See Main c and c lambda, mu = convenient constants relating dSigma and dTau c a,b,c,rhs,soln = temporary workspace used for call to Tridiag c C (BEGIN) start = 1 finish = numTauPts - 1 c initialize the coefficient "bands" in the tridiagonal matrix. c these "bands" remain the same regardless. mu = 0.5*thetanu/dTau lambda = Cnu*dSigma/(dTau**2) do 20010 i=start,finish a(i) = -mu-lambda b(i) = 1.0 + 2.0*lambda c(i) = mu-lambda 20010 continue c Here we will solve the tridiagonal matrix for each "rho" coordinate, c thereby generating solutions for the P values with the same "rho" c but at the next sigma step. First the right-hand-side matrix is c initialized, then the tridiagonal routine is called to solve it, and c solutions are stored back into the P array to record the P values c at the NEXT sigma step. Each "rho" line is assumed to be independent c of the other rho lines, which allows us to solve this as a series of c 2-dimensional finite difference systems (sigma, tau). cvd$ cncall do 20040 j=0,(numRhoPts-1) do 20020 i=start,finish rhs(i) = -mu*P(i-1,j)+P(i,j)+mu*P(i+1,j) 20020 continue call Tridiag(start,finish,a,b,c,rhs,soln,beta,gamma,numTauPts) do 20030 i=start,finish P(i,j) = soln(i) 20030 continue 20040 continue return end c c c SUBROUTINE IBFDdiffract(P, numTauPts, numRhoPts, sigma, + dSigma, dTau, dRho, sumF, + a,b,c,rhs,soln,beta,gamma) c c This routine performs the implicit backward finite difference c method for the diffraction. c implicit none integer numTauPts, numRhoPts real P(0:numTauPts,0:numRhoPts) real sigma, dSigma, dTau, dRho real a(0:numRhoPts),b(0:numRhoPts),c(0:numRhoPts) real rhs(0:numRhoPts),soln(0:numRhoPts) real beta(0:numTauPts),gamma(0:numTauPts) real sumF(0:numRhoPts) real lambda integer i,j, start, finish C VARIABLE DEFINITIONS c See Main. c sumF = summation given by Y.S.L. c lambda = convenient constant relating dSigma, sigma, dRho, dTau c a,b,c...= workspace reserved for Tridiag c (begin) lambda = dTau*dSigma/((1.0+sigma)**2 * dRho**2) c initialize the coefficient "bands" in the tridiagonal matrix. c these "bands" remain the same regardless. a(0) = 0.0 b(0) = 1.0 + lambda/2.0 c(0) = -lambda/2.0 do 30000 j=1,(numRhoPts-2) a(j) = (-1.0 + 1.0/(2.0*j))*lambda/8.0 b(j) = 1.0 + lambda/4.0 c(j) = -(1.0 + 1.0/(2.0*j))*lambda/8.0 30000 continue a(numRhoPts-1) = (-1.0 + 1.0/(2.0*(numRhoPts-1)))*lambda/8.0 b(numRhoPts-1) = 1.0 + lambda/4.0 c(numRhoPts-1) = 0.0 start=0 finish=numRhoPts-1 do 30050 j=start,finish sumF(j) = 0.0 30050 continue c c evaluate the finite differences over "tau line" (lines of constant c tau). Hence the exterior loop is "i" (tau) index. cvd$ cncall do 30250 i=1,(numTauPts-1) c perform the summation step given by Y.S.L. do 30100 j=start,finish sumF(j) = sumF(j) + P(i-1,j) 30100 continue c prep for and solve the system of equations c rhs(0) = P(i,0) - lambda*sumF(0) + lambda*sumF(1) do 30150 j=1,(numRhoPts-2) rhs(j) = P(i,j) + + (1.0 - 1.0/(2.0*j))*lambda*sumF(j-1)/4.0 - + lambda*sumF(j)/2.0 + + (1.0 + 1.0/(2.0*j))*lambda*sumF(j+1)/4.0 30150 continue rhs(numRhoPts-1) = P(i,numRhoPts-1) + + (1.0 -1.0/(2.0*(numRhoPts-1)))*lambda*sumF(numRhoPts-2)/4.0 + - lambda*sumF(numRhoPts-1)/2.0 call Tridiag(start,finish,a,b,c,rhs,soln,beta,gamma,numRhoPts) do 30200 j=start,finish P(i,j) = soln(j) 30200 continue 30250 continue return end c c c SUBROUTINE CNFDabsorb (P, numTauPts, numRhoPts, dSigma, dTau, + Kabsorption, a,b,c,rhs,soln,beta,gamma) c c This routine performs the Crank-Nicolson finite difference method c for the absorption. c c The routine assumes that each element in the array affects c only those elements adjacent to it in retarded time tau, NOT c those adjacent to it radially (in rho). Therefore, each "rho line" c in the array will only affect the corresponding rho line at the c next sigma step, not those adjacent to it. c c The rho lines thus become independent of each other, and allow the c finite differencing to be carried out on each row independently. c Therefore, to conserve memory usage, after computing the solutions c for a particular rho line, the solutions are stored back into the c original P array, taking advantage of the rho line's independence. c implicit none integer numTauPts, numRhoPts real P(0:numTauPts,0:numRhoPts) real dSigma, dTau, Kabsorption real a(0:numTauPts),b(0:numTauPts),c(0:numTauPts) real rhs(0:numTauPts), soln(0:numTauPts) real beta(0:numTauPts),gamma(0:numTauPts) integer start, finish real lambda integer i,j c VARIABLE DEFINITIONS: See Main c and c lambda = a convenient constant relating dSigma and dTau c a,b,c,rhs,soln = temporary workspace used for call to Tridiag c C (BEGIN) start = 1 finish = numTauPts - 1 c initialize the coefficient "bands" in the tridiagonal matrix c these "bands" remain the same regardless lambda = Kabsorption*dSigma/(dTau**2) do 30010 i=start,finish a(i) = -lambda/2.0 b(i) = (1.0 + lambda) c(i) = -lambda/2.0 30010 continue c here we will solve the tridiagonal matrix for each "rho" coordinate, c thereby generating solutions for the P values with the same "rho" c but at the next sigma step. First the right-hand-side matrix is c initialized, then the tridiagonal routine is called to solve it, and c solutions are stored back into the P array to record the P values c at the NEXT sigma step. Each "rho" line is assumed to be independent c of the other rho lines, which allows us to solve this as a series of c 2-dimensional finite differencing systems (sigma, tau). cvd$ cncall do 30040 j=0,(numRhoPts-1) do 30020 i=start,finish rhs(i) = lambda*P(i-1,j)/2.0 + (1-lambda)*P(i,j) + + lambda*P(i+1,j)/2.0 30020 continue call Tridiag(start,finish,a,b,c,rhs,soln,beta,gamma,numTauPts) do 30030 i=start,finish P(i,j) = soln(i) 30030 continue 30040 continue return end c c c SUBROUTINE CNFDrelax (P, numTauPts, numRhoPts, dSigma, dTau, + Cnu, thetanu, a,b,c,rhs,soln,beta,gamma) c c This routine performs the Crank-Nicolson finite difference method c for a single relaxation process. c c The routine assumes that each element in the array affects c only those elements adjacent to it in retarded time tau, NOT c those adjacent to it radially (in rho). Therefore, each "rho line" c in the array will only affect the corresponding rho line at the c next sigma step, not those adjacent to it. c c The rho lines thus become independent of each other, and allow the c finite differencing to be carried out on each row independently. c Therefore, to conserve memory usage, after computing the solutions c for a particular rho line, the solutions are stored back into the c original P array, taking advantage of the rho line's independence. c implicit none integer numTauPts, numRhoPts real P(0:numTauPts,0:numRhoPts) real dSigma, dTau, Cnu, thetanu real a(0:numTauPts),b(0:numTauPts),c(0:numTauPts) real rhs(0:numTauPts), soln(0:numTauPts) real beta(0:numTauPts),gamma(0:numTauPts) integer start, finish real lambda, mu integer i,j c VARIABLE DEFINITIONS: See Main c and c lambda, mu = convenient constants relating dSigma and dTau c a,b,c,rhs,soln = temporary workspace used for call to Tridiag c C (BEGIN) start = 1 finish = numTauPts - 1 c initialize the coefficient "bands" in the tridiagonal matrix c these "bands" remain the same regardless mu = 0.5*thetanu/dTau lambda = 0.5*Cnu*dSigma/(dTau**2) do 30010 i=start,finish a(i) = -mu-lambda b(i) = (1.0 + 2.0*lambda) c(i) = mu-lambda 30010 continue c here we will solve the tridiagonal matrix for each "rho" coordinate, c thereby generating solutions for the P values with the same "rho" c but at the next sigma step. First the right-hand-side matrix is c initialized, then the tridiagonal routine is called to solve it, and c solutions are stored back into the P array to record the P values c at the NEXT sigma step. Each "rho" line is assumed to be independent c of the other rho lines, which allows us to solve this as a series of c 2-dimensional finite differencing systems (sigma, tau). cvd$ cncall do 30040 j=0,(numRhoPts-1) do 30020 i=start,finish rhs(i) = -c(1)*P(i-1,j) + (1-2.0*lambda)*P(i,j) + - a(1)*P(i+1,j) 30020 continue call Tridiag(start,finish,a,b,c,rhs,soln,beta,gamma,numTauPts) do 30030 i=start,finish P(i,j) = soln(i) 30030 continue 30040 continue return end c c c SUBROUTINE CNFDdiffract(P, numTauPts, numRhoPts, sigma, + dSigma, dTau, dRho, sumC, + a,b,c,rhs,soln,beta,gamma) implicit none integer numTauPts, numRhoPts real P(0:numTauPts,0:numRhoPts) real sigma, dSigma, dTau, dRho real a(0:numRhoPts),b(0:numRhoPts),c(0:numRhoPts) real rhs(0:numRhoPts),soln(0:numRhoPts) real beta(0:numTauPts),gamma(0:numTauPts) real sumC(0:numRhoPts) real lambda integer i,j, start, finish c c This routine performs the CRANK-NICOLSON finite difference method C for the diffraction. c C VARIABLE DEFINITIONS c See Main c sumC = summation given by Yang-Sub Lee c lambda = convenient constant relating dSigma, sigma, dRho, dTau c a,b,c...= workspace for Tridiag c (begin) lambda = dTau*dSigma/((1.0+sigma-dSigma/2.0)**2 * dRho**2) c initialize the coefficient "bands" in the tridiagonal matrix. c these "bands" remain the same regardless. a(0) = 0.0 b(0) = 1.0 + lambda/4.0 c(0) = -lambda/4.0 do 40000 j=1,(numRhoPts-2) a(j) = -(1.0 - 1.0/(2.0*j))*lambda/16.0 b(j) = 1.0 + lambda/8.0 c(j) = -(1.0 + 1.0/(2.0*j))*lambda/16.0 40000 continue a(numRhoPts-1) = -(1.0 - 1.0/(2.0*(numRhoPts-1)))*lambda/16.0 b(numRhoPts-1) = 1.0 + lambda/8.0 c(numRhoPts-1) = 0.0 start=0 finish=numRhoPts-1 do 40050 j=start,finish sumC(j) = 0.0 soln(j) = 0.0 40050 continue c c Evaluate the finite differences over "tau line" (lines of c consitant tau). Hence the exterior loop is "i" (tau) index. cvd$ cncall do 40250 i=1,(numTauPts-1) c perform the summation step given by Y.S.L. do 40100 j=start,finish sumC(j) = sumC(j) + soln(j) 40100 continue c prep for and solve the system of equations c rhs(0) = 2.0*P(i,0) - lambda*sumC(0)/2.0 + lambda*sumC(1)/2.0 do 40150 j=1,(numRhoPts-2) rhs(j) = 2.0*P(i,j) + + (1.0 - 1.0/(2.0*j))*lambda*sumC(j-1)/8.0 - + lambda*sumC(j)/4.0 + + (1.0 + 1.0/(2.0*j))*lambda*sumC(j+1)/8.0 40150 continue rhs(numRhoPts-1) = 2.0*P(i,numRhoPts-1) + + (1.0 -1.0/(2.0*(numRhoPts-1)))*lambda*sumC(numRhoPts-2)/8.0 + - lambda*sumC(numRhoPts-1)/4.0 call Tridiag(start,finish,a,b,c,rhs,soln,beta,gamma,numRhoPts) do 40200 j=start,finish P(i,j) = soln(j) - P(i,j) 40200 continue 40250 continue return end c c c SUBROUTINE Nonlinear (P, sigma, dSigma, dTau, tau, t, Pold, + Knonlinearity, numTauPts, numRhoPts, + outfile, dSigShock) c c This routine adds the nonlinear effects to the P array c on a rho-by-rho basis -- that is, it computes the nonlinear effects c for each "rho line" in the P array independently. c c the "t" array is the distorted time base c implicit none real sigma,dSigma, dTau, Knonlinearity integer numTauPts,numRhoPts real P(0:numTauPts,0:numRhoPts) real Pold(0:numTauPts) real tau(0:numTauPts), t(0:numTauPts) real disto, dpdt, dpmax, dSigShock integer i,j,k, outfile C VARIABLE DECLARATIONS c see Main c t = distorted time base c Pold = old pressure waveform for a given rho-line c disto = distortion factor for nonlinear effects c c (begin) disto=Knonlinearity*ALOG(1.0+dSigma/(1.0+sigma-dSigma)) c initialize the ends of t array. The rest will be automatically done c during the continuity test below. These ends are never modified. t(0) = tau(0) t(numTauPts) = tau(numTauPts) dpmax=0.0 do 61000 j=0,numRhoPts-1 c distort the phase of the wavelet to check that continuity c is maintained. Pold(0)=P(0,j) do 60020 i=1,numTauPts Pold(i)=P(i,j) t(i) = tau(i) - Pold(i)*disto if (t(i).LE.t(i-1)) go to 62000 60020 continue c Resample the distorted waveform with linear interpolation c onto the uniform time base k=1 do 60070 i=1,(numTauPts-1) do 60060 while (tau(i).GT.t(k)) k=k+1 if (k.GT.numTauPts) then do 60030 k=i,numTauPts P(k,j)=P(numTauPts,j) 60030 continue go to 61000 endif 60060 continue dpdt=(Pold(k)-Pold(k-1))/(t(k)-t(k-1)) if (dpdt.GT.dpmax) dpmax=dpdt P(i,j)=Pold(k-1)+dpdt*(tau(i)-t(k-1)) 60070 continue 61000 continue dSigShock=(1.0+sigma)/(Knonlinearity*dpmax) return c The discontinuity error message will halt the program 62000 write(*,*) 'Discontinuity: at sigma=',sigma write(*,*) ' smaller dSigma step required' write(outfile,*) 'Discontinuity: at sigma=',sigma stop end c c c SUBROUTINE Tridiag(start, finish, a,b,c, rhs, + soln, beta, gamma, n) c This routine solves the tridiagonal system of (FINISH-START+1) c equations. The arrays A,B,C are the diagonals of the matrix c defining the tridiagonal system. START and FINISH are indices c for the beginning and ending corners that define a SUBMATRIX c within the original matrix. The routine solves the submatrix c using the THOMAS ALGORITHM and stores the results into the c SOLN array. c c c BETA and GAMMA are used as temporary workspace. c implicit none integer start, finish, n, i, j real a(0:n), b(0:n), c(0:n), rhs(0:n), soln(0:n) real beta(0:n), gamma(0:n) beta(start) = b(start) gamma(start)= rhs(start)/beta(start) do 90000 i=(start+1),finish beta(i) = b(i) - a(i)*c(i-1)/beta(i-1) gamma(i) = (rhs(i) - a(i)*gamma(i-1))/beta(i) 90000 continue soln(finish) = gamma(finish) do 90020 i= (finish-1),start,-1 soln(i) = gamma(i) - c(i)*soln(i+1)/beta(i) 90020 continue return end