* 789012345678901234567890123456789012345678901234567890123456789012 * 1 2 3 4 5 6 7 program KeplerAssim * ================================================================== * Assimillation of planetary observations to determine orbital * parameters. Fortran 77. * Observations and other data are read-in from a file. A run of * this code does one of four things: * 1. XSEC plot cross sections through the cost function. * 2. ANAL analyse observations for their best-fit orbital params. * 3. RECO Turn alt/azi observations in input file into RA/dec * observations in a new input file (for trials only). * 4. TRAJ Use params in input file to generate planetary * trajectories. * See documentation "www.met.rdg.ac.uk/~ross/Astronomy/Kepler.html" * for further info. * ================================================================== * Development history: * Ross Bannister, February 2001 to February 2003 * Include modifications as recommended by reviewers of paper * 1. Option for second derivs of forward model in Hessian calc * (direct inversion maintained), December 2002. * 2. Use double precision numbers, February 2003. * 3. Include Gaussian ellimination, February 2003 * ================================================================== * Declare variables * ================= implicit none * ----- The No of observations (#alt+#azi+#RA+#dec measurements) is * Nobs. This should be <= maxobs integer maxobs,Nobs parameter (maxobs=50) * ----- Variables that hold information about the obs. integer tpe(maxobs) real*8 y0(maxobs),Error(maxobs),Time(maxobs),long,lat,x(6) real*8 y(maxobs) character*8 filen character*4 rmode character*12 fname logical verbose,CorrectRefract,FinalRun,d2Hess * ----- Variables that are used in the analysis integer i,j,k,iter,maxiter,partial real*8 pi,ms,crit,CostFunc,amax,amin,imax,imin,range(6),x0(6) real*8 H(maxobs,6),Hess(6,6),Hess1(6,6),Hess2(6,6),Hd2(maxobs,6,6) real*8 InvH(6,6),Grad(6) real*8 Temp1(maxobs),Temp2(maxobs),Temp3(maxobs),Temp4(maxobs) real*8 Lambda(6,6),U(6,6),deltay(maxobs),deltax(6) real*8 total,low,high,v(6),trace,CostJ,son,lim * ----- Variables associated with the forward model integer MonthDays(0:12),zeroYr real*8 RefTime,NumDays,obl,Epoch,aE,eE,iE,Omega1E,omega2E,LE * Initialze Data * ============== data MonthDays/0,31,59,90,120,151,181,212,243,273,304,334,365/ pi=4.0*atan(1.0) * Zero of time will be 1st Jan 1990 zeroYr=1990 * Reference time for alt\azi conversion RefTime=NumDays(zeroYr,1990,3,21,0,0,MonthDays) obl=pi*(23.0+(26.0+21.0/60.0)/60.0)/180.0 * Epoch is called t0 in many subroutines Epoch=NumDays(zeroYr,2000,1,1,0,0,MonthDays) * Earth parameters data aE,eE,iE,Omega1E,omega2E,LE/1.00000011, 0.01671022, 0.00005, : -11.26064, 102.94719, 100.46435/ * Read-in observational data, incl errs and long \ lat, etc. * ========================================================== call ReadObs (tpe,y0,Error,Time,Nobs,maxobs,zeroYr,long,lat,x, : rmode,verbose,amin,amax,imin,imax,crit,maxiter, : partial,CorrectRefract,d2Hess,MonthDays,fname) * y0 array contains the observations * Error contains the associated errors * Time contains the time of observations * x is the inital guess at the parameters * rmode is the run mode * amin,amax is the range of parameter a (similarly for imin,imax) if (verbose) then print*,'=======================================================' print*,'Astronomical Data Assimmilation' print*,'Ross Bannister, February/April 2001 & July/August 2002' print*,'University of Reading, UK' print*,'=======================================================' print* print*,'longitude =',long print*,'latitude =',lat print*,'Number of observation pairs =',Nobs/2 print*,'Initial guess: ',(x(k),k=1,6) print*,'Range of parameter a: ',amin,amax print*,'Range of parameter i: ',imin,imax print*,'Run mode = ',rmode endif * 'IF' STATEMENTS THAT EXECUTE THE CHOSEN RUN MODE OF THE CODE * ============================================================ * ------------------------------------------------------------------ if (rmode.eq.'XSEC') then * ------------------------- if (verbose) print*,'RUN MODE: Compute cost fn cross sections' * Probe some of physical parameter space to find out how J behaves * Define limits of the varibles about the x already defined above data range/5.0,0.4,180.,180.,180.,180./ do k=1,6 x0(k)=x(k) enddo do k=1,6 * Adjust parameter k in turn * Open a file for the output write (filen,'(A5,I1)') 'Param',k if (verbose) print*,filen open (12,file=filen) do i=0,150 * Set range for parameter to vary over if (k.eq.1) then x0(k)=amin+real(i)*(amax-amin)/150. elseif (k.eq.3) then x0(k)=imin+real(i)*(imax-imin)/150. else x0(k)=x(k)+range(k)*real(i-75)/75. endif call CheckRange (x,verbose) * Find the projected observations with this set of parameters call AstroFM (x0,aE,eE,iE,Omega1E,omega2E,LE,Time,Epoch, : long,lat,RefTime,tpe,Nobs,maxobs,y,obl) * Predicted observations are in y * Find the cost function write (12,202) x0(k),CostFunc(y,y0,Error,Nobs,maxobs) enddo close (12) x0(k)=x(k) enddo elseif (rmode.eq.'ANAL') then * ----------------------------- * ------------------------------------------------------------------ * FinalRun allows a subset of the code within the ANAL loop to be * rerun after the analysis has been computed (to help compute * analysis errors). Set to true when converged FinalRun=.false. if (verbose) print*,'RUN MODE: Perform least squares minimization' * Open a file for the output open (12,file='Convergence') write (12,203) '#============================================='// : '===================================' write (12,203) '# No. a e i Omega '// : 'omega-bar Epsilon Inc. Cost' write (12,203) '#---------------------------------------------'// : '-----------------------------------' iter=0 * Output the initial parameters and cost function call AstroFM (x,aE,eE,iE,Omega1E,omega2E,LE,Time,Epoch, : long,lat,RefTime,tpe,Nobs,maxobs,y,obl) write (12,201) iter,(x(i),i=1,6),ms, : CostFunc(y,y0,Error,Nobs,maxobs) 10 continue * This is the start of the minimization loop * ========================================== if (verbose) print*,'--------------------------------------------' * (x is the linearization point) if (FinalRun) then print*,'##### PERFORMING ERROR ANALYSIS #####' else iter=iter+1 print*,'##### STARTING ITERATION NO. ',iter,' #####' endif 20 continue if (.not.FinalRun) then if (iter.le.partial) then * Do a partial minimization over inclination (parameter 3) call PartialMin (3,imin,imax,x,aE,eE,iE,Omega1E,omega2E,LE, : Time,Epoch,long,lat,RefTime,tpe,Nobs,y,obl,y0,Error, : maxobs,verbose) * Do a partial minimization over semi-major axis (parameter 1) call PartialMin (1,amin,amax,x,aE,eE,iE,Omega1E,omega2E,LE, : Time,Epoch,long,lat,RefTime,tpe,Nobs,y,obl,y0,Error, : maxobs,verbose) * Put the parameters into the appropriate range call CheckRange (x,verbose) endif * Calculate the projected obs from the current estimate of x * Use the forward model if (verbose) print*,'Calculating the forward model' call AstroFM (x,aE,eE,iE,Omega1E,omega2E,LE,Time,Epoch, : long,lat,RefTime,tpe,Nobs,maxobs,y,obl) * Array y contains projected observations (y0 contains actual obs) if (verbose) then print*,'Here are the observations and projected observations' do k=1,Nobs print*,int(real(k+1)/2.),180.*y0(k)/pi,180.*y(k)/pi enddo endif * Find deltay do k=1,Nobs deltay(k)=y0(k)-y(k) enddo * Due to the periodic nature of angles, deltay calculated in this * way may be misleading. Correct for this here. do k=1,Nobs if (tpe(k).eq.1) then * The measurement is of RA or dec (within range 0-2pi) * Need to possibly adjust RA if (mod(k,2).eq.1) call PutInRange (deltay(k),pi,0) else * The measurement is of alt or azi (within range 0-2pi) * Need to possibly adjust azi if (mod(k,2).eq.1) call PutInRange (deltay(k),pi,0) endif enddo endif * Linearize the forward model (for the Jacobean) * This is where we need real*8 instead of real variables if (verbose) print*,'Linearizing forward model' call LinearFM (x,aE,eE,iE,Omega1E,omega2E,LE,Time,Epoch,long,lat, : RefTime,tpe,Nobs,maxobs,H,Temp1,Temp2,obl) if (verbose) then * Output the forward model print*,'----- Jacobean ------------------------------------' do k=1,Nobs print 204,k,(H(k,i),i=1,6) enddo print*,'---------------------------------------------------' endif * -----Calculate the Hessian in parameter space-------------------- if (verbose) print*,'Calculating Hessian - 1st derivs of H' * Contribution from the 1st derivatives of H (H(T)R(-1)H) do i=1,6 do j=1,6 if (j.ge.i) then * Upper triangular and diagonal terms total=0. do k=1,Nobs total=total+H(k,i)*H(k,j)/Error(k) enddo Hess1(i,j)=total else * Lower diagonal terms Hess1(i,j)=Hess1(j,i) endif Hess(i,j)=Hess1(i,j) enddo if (verbose) print 200,(Hess1(i,j),j=1,6) enddo if (d2Hess) then * Contribution from the 2nd derivatives of H (if required) if (verbose) print*,'Calculating Hessian - 2nd derivs of H' call SecondH (x,aE,eE,iE,Omega1E,omega2E,LE,Time,Epoch,long,lat, : RefTime,tpe,Nobs,maxobs,Hd2,obl,Temp1,Temp2,Temp3, : Temp4) * Compute the rest of the Hessian do i=1,6 do j=1,6 if (j.ge.i) then * Upper triangular and diagonal terms total=0. do k=1,Nobs total=total-Hd2(k,i,j)*deltay(k)/Error(k) enddo Hess2(i,j)=total else * Lower diagonal terms Hess2(i,j)=Hess2(j,i) endif Hess(i,j)=Hess(i,j)+Hess2(i,j) enddo if (verbose) print 200,(Hess2(i,j),j=1,6) enddo if (verbose) then * Display full Hessian print*,'Total Hessian' do i=1,6 print 200,(Hess(i,j),j=1,6) enddo endif endif if (.not.FinalRun) then * -----Calculate minus the gradient in parameter space------------ do i=1,6 total=0. do k=1,Nobs total=total+H(k,i)*deltay(k)/Error(k) enddo Grad(i)=total enddo * -----Solve equations to find increment in x-space--------------- * Solve: Hessian deltax = -gradient * Use Gaussian ellimination - the increment is deltax * Hess1 and v in this call are now just workspace variables call GaussEll (Hess,Hess1,Grad,v,deltax,6,6) * -----Calculate updated x (and find increment size)-------------- ms=0. do i=1,6 ms=ms+deltax(i)*deltax(i) * Find the updated value of x x(i)=x(i)+deltax(i) enddo * Put the parameters into the appropriate range call CheckRange (x,verbose) CostJ=CostFunc(y,y0,Error,Nobs,maxobs) write (12,201) iter,(x(i),i=1,6),ms,CostJ print 201,iter,(x(i),i=1,6),ms,CostJ * Check that a and i remain in the specified ranges if ((x(1).lt.amin).or.(x(1).gt.amax)) then print*,'Fatal error from main program' print*,'param 1 is not in the specified range' print*,'value : ',x(1) print*,'Range : ',amin,amax stop endif if ((x(3).lt.imin).or.(x(3).gt.imax)) then print*,'Fatal error from main program' print*,'param 3 is not in the specified range' print*,'value : ',x(3) print*,'Range : ',imin,imax stop endif endif * -----Have we converged?------------------------------------------- if ((.not.FinalRun).and.(ms.gt.crit).and.(iter.lt.maxiter)) : goto 10 if ((ms.gt.crit).and.verbose) : print*,'Premature exit - maximum No. of iterations reached' * Have completed the iterations - Set FinalRun flag for last run * through (need to repeat some of algorithm - for error analysis) if (.not.FinalRun) then FinalRun=.true. goto 10 endif * End of iterations write (12,203) '#---------------------------------------------'// : '-----------------------------------' * -----Find the inverse Hessian (for error calculation)------------- * A final and partial trip round the loop (FinalRun true) above * calculated the Hessian only if (verbose) print*,'Diagonalizing' lim=0.00001 call Jacobi (Hess,Lambda,U,6,6,lim,trace) * The eigenvalues are the diagonal elements of Lambda * The eigenvectors are the rows of U if (verbose) then print*,'Here are the eigenvalues of the Hessian' print 200,(Lambda(i,i),i=1,6) print*,'Here is the matrix of eigenvectors' do i=1,6 print 200,(U(i,j),j=1,6) enddo endif * Check for eigenvalues that show Hessian is not pos. definite do i=1,6 if (Lambda(i,i).le.0.) : print*,'Warning - bad eigenvalue No. ',i,' (',Lambda(i,i),')' enddo * Compute the inverse Hessian matrix, U(T)lambda(-1)U do i=1,6 do j=1,6 total=0. do k=1,6 if (Lambda(k,k).ne.(0.)) : total=total+U(k,i)*U(k,j)/Lambda(k,k) enddo InvH(i,j)=total enddo enddo print*,'---------------------------------------------------------' write (12,203) '# Analysis errors (negative numbers have negative' : //' curvatures)' write (12,205) '# ',(son(InvH(i,i))*sqrt(abs(InvH(i,i))),i=1,6) write (12,203) '#============================================='// : '===================================' write (12,203) '# Hessian eigenvalues' write (12,206) '# ',(Lambda(i,i),i=1,6) write (12,203) '#============================================='// : '===================================' close (12) print*,'Convergence and analysis details are in file Convergence' elseif (rmode.eq.'RECO') then * ----------------------------- if (verbose) then print*,'RUN MODE: Reconfigure input file' print*,'Any alt/azi details will be transformed to RA/dec' endif call Reconfig (fname,zeroYr,RefTime,MonthDays) elseif (rmode.eq.'TRAJ') then * ----------------------------- if (verbose) then print*,'RUN MODE: Generate RA/dec trajectories' print*,'Any alt/azi details will be transformed to RA/dec' endif * Find the min and max time call TimeRange (Time,Nobs,low,high,maxobs) * Open two files 12: model trajectory 'low' to 'high' (TrajCont) * 13: obs and model at obs times (TrajPoints) open (12,file='TrajCont') open (13,file='TrajPoints') * Deal with the continuous trajectory (250 points) if (verbose) print*,'Continuous trajectory calculation' call TrajCont (low,high,250,x,aE,eE,iE,Omega1E,omega2E,LE,Epoch, : obl,12) * Deal with the obs and projected obs if (verbose) print*,'Obs and predicted obs calculation' call TrajPoint (Time,tpe,y0,Nobs,maxobs,reftime,long,lat,x,aE,eE, : iE,Omega1E,omega2E,LE,Epoch,obl,13) close (12) close (13) if (verbose) print*,'Output in files TrajCont and TrajPoints' endif * ----- print*,'Program finished' 200 format (6F14.5) 201 format (I4,F8.3,F9.4,5F10.4,F10.1) 202 format (F12.6,F16.2) 203 format ((A)) 204 format (I6,F10.2,F10.3,F14.7,F15.8,F17.10,F14.7) 205 format ((A),F8.3,F9.4,4F10.4) 206 format ((A),6F13.3) end * ================================================================== * ================================================================== real*8 function CostFunc (y,y0,Error,Nobs,maxobs) * This routine will actually calculate the value of the cost * function (x-representation, non linear) * y0 observations * y projected observations implicit none integer maxobs,i,Nobs real*8 y(maxobs),y0(maxobs),Error(maxobs),total,diff,pi pi=4.*atan(1.) total=0. do i=1,Nobs diff=abs(y0(i)-y(i)) * Put into correct range call PutInRange (diff,pi,0) total=total+(diff)**2./Error(i) enddo CostFunc=total/2. return end * ================================================================== subroutine ReadObs (tpe,RawObs,Error,Time,c,maxobs,zeroYr,long, : lat,x,rmode,verbose,amin,amax,imin,imax,crit, : maxiter,partial,CorrectRefract,d2Hess, : MonthDays,fname) * This subroutine will read a AstroObs file * All angles output are in radians (except long and lat - degrees) implicit none integer maxobs,tpe(maxobs),c,c1,t,d,Yr,Mo,Dy,Hr,Mn,zeroYr,s,l integer MonthDays(0:12),switch,maxiter,partial real*8 RawObs(maxobs),Time(maxobs),Error(maxobs),err,NumDays,long real*8 lat,pi,deg,x(6),amin,amax,imin,imax,crit character*75 temp character*12 fname character*4 rmode logical verbose,CorrectRefract,d2Hess pi=4.*atan(1.) print*,'Please enter the input filename (e.g. AstroObs)' read*,fname print*,fname open (12,file=fname) 12 read (12,100,err=13) temp if (index(temp,'END OF HEADER').eq.0) goto 12 goto 14 13 print*,'Error: No end to input file header.' stop 14 read (12,100,err=13) temp read (12,100,err=13) temp * Longitude read (12,*,end=11) s,d,Mn long=real(s)*(real(d)+real(Mn)/60.) read (12,100) temp * Latitude read (12,*,end=11) s,d,Mn lat=real(s)*(real(d)+real(Mn)/60.) read (12,100) temp * First guess at solution read (12,*) (x(l),l=1,6) read (12,100) temp * Run mode (XSEC, ANAL, RECO, TRAJ) read (12,100) rmode read (12,100) temp * Verbose read (12,*) switch verbose=(switch.ne.0) read (12,100) temp * Ranges (length of semi-major axis) read (12,*) amin,amax read (12,100) temp * (inclination) read (12,*) imin,imax read (12,100) temp * Check that the first guess is within this range if ((x(1).lt.amin).or.(x(1).gt.amax)) then print*,'Fatal error from subroutine ReadObs' print*,'Initial guess for param 1 is not in the specified range' print*,'Initial guess : ',x(1) print*,'Range : ',amin,amax stop endif if ((x(3).lt.imin).or.(x(3).gt.imax)) then print*,'Fatal error from subroutine ReadObs' print*,'Initial guess for param 3 is not in the specified range' print*,'Initial guess : ',x(3) print*,'Range : ',imin,imax stop endif * Convergence cutoff in length of increment read (12,*) crit read (12,100) temp * Maximum number of iterations read (12,*) maxiter read (12,100) temp * Number of iterations performing partial minimization read (12,*) partial read (12,100) temp * Refraction correction read (12,*) switch CorrectRefract=(switch.ne.0) read (12,100) temp * Include 2nd derriv of forward model in Hessian? read (12,*) switch d2Hess=(switch.ne.0) read (12,100) temp c=-1 10 read (12,100,end=11) temp read (12,*,end=11) t * Increment the counter c=c+2 c1=c+1 if (c1.gt.maxobs) then print*,'TOO MANY OBSERVATIONS' print*,'INCREASE maxobs' stop endif * Store the type flag tpe(c)=t tpe(c1)=t * Find the time read (12,*,end=11) Yr,Mo,Dy,Hr,Mn Time(c)=NumDays(zeroYr,Yr,Mo,Dy,Hr,Mn,MonthDays) Time(c1)=Time(c) if (t.eq.1) then * The observation is of RA \ dec variety read (12,*,end=11) Hr,Mn,err RawObs(c)=pi*(real(Hr)+real(Mn)/60.)/12. Error(c)=(pi*err/720.)**2. read (12,*,end=11) d,Mn,err RawObs(c1)=pi*(real(d)+real(Mn)/60.)/180. Error(c1)=(pi*err/10800.)**2. else * The observation is of alt \ azi variety read (12,*,end=11) deg,err RawObs(c)=pi*deg/180. Error(c)=(pi*err/180.)**2. read (12,*,end=11) deg,err RawObs(c1)=pi*deg/180. Error(c1)=(pi*err/180.)**2. endif goto 10 11 close (12) c=c+1 100 format ((A)) return end * ================================================================== subroutine ReconFig (fname,zeroYr,RefTime,MonthDays) * This subroutine will reconfigure an AstroObs file * (Change all alt/azi observations to RA/dec type) * All angles output are in radians (except long and lat - degrees) implicit none integer t,d,Yr,Mo,Dy,Hr,Mn,zeroYr,s,l integer MonthDays(0:12),verb,maxiter,partial,switch real*8 Time,err,NumDays,RefTime,long,RA,dec,crit real*8 lat,pi,x(6),amin,amax,imin,imax,alt,azi,erralt,errazi,son character*75 temp character*4 rmode character*12 fname pi=4.*atan(1.) print*,'Opening input file for reading : ',fname open (12,file=fname) print*,'Opening output file for writing : AstroObs1' open (13,file='AstroObs1') do l=1,28 read (12,100) temp write (13,100) temp enddo * Longitude read (12,*,end=11) s,d,Mn long=real(s)*(real(d)+real(Mn)/60.) write (13,*) s,d,Mn read (12,100) temp write (13,100) temp * Latitude read (12,*,end=11) s,d,Mn lat=real(s)*(real(d)+real(Mn)/60.) write (13,*) s,d,Mn read (12,100) temp write (13,100) temp * First guess at solution read (12,*) (x(l),l=1,6) write (13,*) (x(l),l=1,6) read (12,100) temp write (13,100) temp * Run mode (1=cost fn cross sections, 2=minimize cost fn, * 3=reconfigure input file alt/azi -> RA/dec * 4=generate RA/dec trajectory) read (12,100) rmode write (13,100) rmode read (12,100) temp write (13,100) temp * Verbose read (12,*) verb write (13,*) verb read (12,100) temp write (13,100) temp * Ranges (length of semi-major axis) read (12,*) amin,amax write (13,*) amin,amax read (12,100) temp write (13,100) temp * (inclination) read (12,*) imin,imax write (13,*) imin,imax read (12,100) temp write (13,100) temp * Convergence cutoff in length of increment read (12,*) crit write (13,*) crit read (12,100) temp write (13,100) temp * Maximum number of iterations read (12,*) maxiter write (13,*) maxiter read (12,100) temp write (13,100) temp * Number of iterations performing partial minimization read (12,*) partial write (13,*) partial read (12,100) temp write (13,100) temp * Refraction correction read (12,*) switch write (13,*) switch read (12,100) temp write (13,100) temp * Include 2nd derriv of forward model in Hessian? read (12,*) switch write (13,*) switch read (12,100) temp write (13,100) temp 10 read (12,100,end=11) temp write (13,100) temp read (12,*,end=11) t write (13,*) 1 * Find the time read (12,*,end=11) Yr,Mo,Dy,Hr,Mn write (13,*) Yr,Mo,Dy,Hr,Mn Time=NumDays(zeroYr,Yr,Mo,Dy,Hr,Mn,MonthDays) if (t.eq.1) then * The observation is of RA \ dec variety - no need to change read (12,*,end=11) Hr,Mn,err write (13,*) Hr,Mn,err read (12,*,end=11) d,Mn,err write (13,*) d,Mn,err else * The observation is of alt \ azi variety - need to convert read (12,*,end=11) alt,erralt read (12,*,end=11) azi,errazi call localinverse (RA,dec,Time,long,lat,reftime,alt,azi) * NEED TO DEAL WITH THE ERROR err=2. * Deal with the RA s=int(son(RA)) RA=abs(RA) Hr=s*int(RA) RA=60.*(RA-real(abs(Hr))) Mn=int(RA) write (13,*) Hr,Mn,err * Deal with the dec s=int(son(dec)) dec=abs(dec) d=s*int(dec) dec=60.*(dec-real(abs(d))) Mn=int(dec) write (13,*) d,Mn,err endif goto 10 11 close (12) close (13) 100 format ((A)) return end * ================================================================== subroutine LinearFM (x0,aE,eE,iE,Omega1E,omega2E,LE,Time,t0,long, : lat,reftime,tpe,Nobs,maxobs,H,Obs1,Obs2,obl) * Linearize the forward model about the point x0. * The vector x0 has components: * x0(1)=a, x0(2)=e, x0(3)=i, x0(4)=Omega1, x0(5)=omega2, x0(6)=L implicit none integer maxobs,tpe(maxobs),Nobs,param,obs real*8 x0(6),aE,eE,iE,Omega1E,omega2E,LE,Time(maxobs),t0,long,lat real*8 reftime,H(maxobs,6),delta,Obs1(maxobs),Obs2(maxobs),d real*8 obl,x0up(6),x0dn(6) do param=1,6 x0up(param)=x0(param) x0dn(param)=x0(param) enddo delta=0.00005 do param=1,6 * (Vary each parameter in turn) * Need slightly different values for gradient calculation x0up(param)=x0(param)+delta x0dn(param)=x0(param)-delta * For the eccentricity, check the range if (param.eq.2) then if (x0up(param).gt.(1.)) x0up(param)=1. if (x0dn(param).lt.(0.)) x0dn(param)=0. endif d=x0up(param)-x0dn(param) * Call the forward model with a slightly modified set of params call AstroFM (x0up,aE,eE,iE,Omega1E,omega2E,LE,Time,t0, : long,lat,reftime,tpe,Nobs,maxobs,Obs1,obl) call AstroFM (x0dn,aE,eE,iE,Omega1E,omega2E,LE,Time,t0, : long,lat,reftime,tpe,Nobs,maxobs,Obs2,obl) * Restore the working list of parameters x0up(param)=x0(param) x0dn(param)=x0(param) * Calculate the components of the derivative do obs=1,Nobs H(obs,param)=(Obs1(obs)-Obs2(obs))/d enddo enddo return end * ================================================================== subroutine AstroFM (xx,aE,eE,iE,Omega1E,omega2E,LE,Time,t0, : long,lat,reftime,tpe,Nobs,maxobs,ProObs,obl) * Astronomical 'forward model'. * Given a set of orbital parameters, the times, the location on * Earth and the type of observation, this subroutine will generate * the observtions which are expected. implicit none integer maxobs,Nobs,tpe(maxobs),loop,pairs,ind1,ind2 real*8 xx(6),a,e,i,Omega1,omega2,L,aE,eE,iE,Omega1E,omega2E,LE,t0 real*8 Time(maxobs),long,lat,ProObs(maxobs),RA,Dec,reftime real*8 alt,azi,obl,xcel,ycel,zcel,xecl,yecl,zecl * Longitude and latitude should be in degrees * Nobs is actually 2*(number of observation co-ordinates) a=xx(1) e=xx(2) i=xx(3) Omega1=xx(4) omega2=xx(5) L=xx(6) pairs=Nobs/2 do loop=1,pairs ind1=(loop-1)*2+1 ind2=ind1+1 call Orbit (a,e,i,Omega1,omega2,L,aE,eE,iE,Omega1E,omega2E, : LE,Time(ind1),t0,obl,xcel,ycel,zcel,xecl,yecl,zecl) * The used o/p is xcel,ycel,zcel (geocentric equatorial co-ords) if (tpe(ind1).eq.1) then * Observation required is RA \ dec * Calculate RA and dec from xcel,ycel,zcel call Celestial (xcel,ycel,zcel,RA,Dec) ProObs(ind1)=RA ProObs(ind2)=Dec else * Observation required is alt \ azi call local (xcel,ycel,zcel,Time(ind1),long,lat,reftime,alt, : azi) ProObs(ind1)=alt ProObs(ind2)=azi endif enddo return end * ================================================================== subroutine local (xeq,yeq,zeq,time,long,lat,reftime,alt,azi) * Convert to local (alt-azimuth) co-ordinates * Expecting xeq,yeq,zeq (equ. co-ords), long and lat in degrees * alt and azi are output in radians implicit none real*8 xeq,yeq,zeq,time,long,lat,reftime,reflong,radlat real*8 slong,clong,slat,clat,Nhor,Ehor,Zhor,alt,azi,pi,remainder pi=4.0*atan(1.0) * Determine the effective longitude (in radians) reflong=long/360.0+(time-reftime)/0.997269562 reflong=2.0*pi*remainder(reflong) * Workout some trig radlat=pi*lat/180.0 slong=sin(reflong) clong=cos(reflong) slat=sin(radlat) clat=cos(radlat) * Convert to local (horizontal) co-ordinates Nhor=slat*clong*xeq+slat*slong*yeq+clat*zeq Ehor=slong*xeq-clong*yeq Zhor=-clat*clong*xeq-clat*slong*yeq+slat*zeq * Calculate the altitude (radians) alt=atan(Zhor/sqrt(Ehor*Ehor+Nhor*Nhor)) * Calculate the azimuth (radians) if (Nhor.eq.0) then if (Ehor.gt.0) then azi=pi/2. else azi=3.*pi/2. endif else azi=atan(Ehor/Nhor) if (Nhor.lt.(0.0)) azi=azi-pi endif if (azi.lt.0) azi=azi+2.*pi return end * ================================================================== subroutine localinverse (RA,dec,time,long,lat,reftime,alt,azi) * Convert from local horizontal co-ords to geocentric equatorial * celestial co-ords. * Main inputs : alt and azi in degrees * Main outputs : RA in hours and dec in degrees * Other intputs : time,long,lat,reftime implicit none real*8 xeq,yeq,zeq,time,long,lat,reftime,reflong,radlat real*8 slong,clong,slat,clat,Nhor,Ehor,Zhor,alt,azi,pi,remainder real*8 RA,dec,altrad,azirad pi=4.0*atan(1.0) * Convert alt and azi to radians altrad=pi*alt/180. azirad=pi*azi/180. * Calculate the local Cartesian co-ordinates from alt and azi Nhor=cos(altrad)*cos(azirad) Ehor=cos(altrad)*sin(azirad) Zhor=sin(altrad) * Determine the effective longitude (in radians) reflong=long/360.0+(time-reftime)/0.997269562 reflong=2.0*pi*remainder(reflong) * Workout some trig radlat=pi*lat/180.0 slong=sin(reflong) clong=cos(reflong) slat=sin(radlat) clat=cos(radlat) * Convert from local (horizontal) co-ordinates to equatorial co-ords xeq=slat*clong*Nhor+slong*Ehor-clat*clong*Zhor yeq=slat*slong*Nhor-clong*Ehor-clat*slong*Zhor zeq=clat*Nhor+slat*Zhor * Calculate RA and dec from the equatorial co-ordinates call Celestial (xeq,yeq,zeq,RA,dec) * Convert RA into hours and dec into degrees RA=12.*RA/pi dec=180.*dec/pi return end * ================================================================== subroutine Orbit (a,e,i,Omega1,omega2,L,aE,eE,iE,Omega1E,omega2E, : LE,t,t0,obl,xcel,ycel,zcel,xecl,yecl,zecl) * Returns the xcel,ycel,zcel (geocentric celestial co-ordinates (can * compute the RA/dec from these)) and xecl,yecl,zecl (ecliptic * co-ordinates) given the orbital parameters of a planet and the * Earth at the chosen time. * Ross Bannister, December 2000 * Input : a,e,i,Omega1,omega2,L (orbital parameters of planet) * t,t0 (time and Epoch in days) * aE,eE,iE,Omega1E,omega2E,LE (orbital parameters of Earth) * Output : Geocentric celestial co-ordinates (xcel,ycel,zcel) * Ecliptic co-ordinates (xecl,yecl,zecl) implicit none real*8 a,e,i,Omega1,omega2,L,t,t0,aE,eE,iE,Omega1E,omega2E,LE real*8 xp,yp,zp,xE,yE,zE,pi,pi2,mymod,Ldot,LEdot,obl real*8 Lrad,Omega1rad,omega2rad,irad,LErad,Omega1Erad,omega2Erad real*8 iErad,con,xcel,ycel,zcel,xecl,yecl,zecl pi=4.0*atan(1.) pi2=2.*pi * con calculates Ldot from a (Ldot=con/a**1.5) con=129597762. Ldot=con/a**1.5 LEdot=con/aE**1.5 * Convert some parameters to radians Lrad=pi*(L+Ldot*(t-t0)/131490000.0)/180.0 Lrad=mymod(Lrad,pi2) Omega1rad=pi*Omega1/180.0 omega2rad=pi*omega2/180.0 irad=pi*i/180.0 LErad=pi*(LE+LEdot*(t-t0)/131490000.0)/180.0 LErad=mymod(LErad,pi2) Omega1Erad=pi*Omega1E/180.0 omega2Erad=pi*omega2E/180.0 iErad=pi*iE/180.0 * Position of planet xp=0 yp=0 zp=0 * Modified planetary co-ordinates call PosPlanetary (a,omega2rad,Lrad,e,Omega1rad,xp,yp) * Convert to ecliptic co-ordinates call ecliptic (xp,yp,zp,Omega1rad,irad) xecl=xp yecl=yp zecl=zp * Position of Earth xE=0 yE=0 zE=0 * Modified planetary co-ordinates call PosPlanetary (aE,omega2Erad,LErad,eE,Omega1Erad,xE,yE) * Convert to ecliptic co-ordinates call ecliptic (xE,yE,zE,Omega1Erad,iErad) * Geocentric co-ordinates of the planet call geo (xp,yp,zp,xE,yE,zE) * Geocentric equatorial co-ordinates call geoequ (yp,zp,obl) * xcel,ycel, and zcel are output xcel=xp ycel=yp zcel=zp return end * ================================================================== subroutine PosPlanetary (a,longperi,mlong,e,longan,x,y) * To calculate the position of a planet in the modified planetary * co-ordinates implicit none integer count real*8 a,longperi,mlong,e,longan real*8 x,y,M,E1,E2,error,value real*8 deriv,change,w,r,argperi,mymod,pi,pi2,Kepler real*8 Kderiv * Determine pi pi=4.0*atan(1.0) pi2=2.*pi * Calculate the mean anomaly, eq. (7) M=mymod(mlong-longperi,pi2) * Use Newton-Raphson procedure to calculate the eccentric anomaly error=0.00001 * First guess at answer E2=M count=0 10 E1=E2 count=count+1 value=Kepler(M,E1,e) deriv=Kderiv(E1,e) change=value/deriv * Estimate better value do eccentric anomaly E2=E1-change if (.not.(ABS(change).lt.error).or.(count.gt.50)) goto 10 * Calculate the true anomaly w=2*atan(sqrt((1+e)/(1-e))*tan(E2/2)) * Distance from occupied focus (AU) r=a*(1-e*cos(E2)) * Calculate the argument of the perihelion argperi=longperi-longan * Evaluate x,y in the modified planetary co-ordinate system (eq. 8) x=r*cos(w+argperi) y=r*sin(w+argperi) return end * ================================================================== real*8 function Kepler(M,E,ecc) * Kepler's equation implicit none real*8 M,E,ecc Kepler=M-E+ecc*sin(E) return end real*8 function Kderiv(E,ecc) * Derivative of Kepler's equation implicit none real*8 E,ecc Kderiv=ecc*COS(E)-1 return end * ================================================================== subroutine ecliptic(x,y,z,longan,i) * Convert from 'modified planetary' to 'ecliptic' co-ordinates implicit none real*8 x,y,z,xx,yy,zz,longan,i,csl real*8 snl,csi,sni * Work out the longitude of the ascending node and the inclination * Calculate some trig functions csl=cos(longan) snl=sin(longan) csi=cos(i) sni=sin(i) * Do the transformation (z=0, so it is not included in the formula) xx=csl*x-snl*csi*y yy=snl*x+csl*csi*y zz=sni*y x=xx y=yy z=zz return end subroutine geo(x,y,z,Earthx,Earthy,Earthz) * Ecliptic to geocentric co-ordinates implicit none real*8 x,y,z,Earthx,Earthy,Earthz x=x-Earthx y=y-Earthy z=z-Earthz return end subroutine geoequ(y,z,obl) * Geocentric to equatorial geocentric co-ordinates * (The x-co-ordinate is not affected) implicit none real*8 y,z,obl,cs,sn,yy,zz cs=cos(obl) sn=sin(obl) yy=cs*y-sn*z zz=sn*y+cs*z y=yy z=zz return end * ================================================================== subroutine Celestial (x,y,z,RA,dec) * IN : x,y,z of celestial object in geocentric equatorial co-ords * OUT: RA (in rads) and dec (in rads) implicit none real*8 x,y,z,RA,dec,pi pi=4.0*atan(1.0) * Calculate the RA in radians if (x.eq.0) then if (y.gt.0) then RA=pi/2.0 else RA=3.0*pi/2.0 endif else RA=atan(y/x) if (x.lt.0) RA=RA+pi if (RA.lt.0) RA=RA+2.0*pi endif * Calculate the dec in radians dec=atan(z/sqrt(x*x+y*y)) return end * ================================================================== * ===== Routines associated with date conversion =================== real*8 function NumDays (zeroYr,Yr,Mo,Dy,Hr,Mi,MonthDays) * Function to calculate the number of days since the zero of time implicit none integer zeroYr,Yr,Mo,Dy,Hr,Mi,MonthDays(0:12),year real*8 days logical leapyr * Check that the time of interest is > zero time if (Yr.lt.zeroYr) then print 101,'YOU MUST CHOOSE A TIME WHICH IS AT OR AFTER',zeroYr print 101,'USING TIME 00:00 GMT Jan 01 ',zeroYr days=0.0 else * Consider the number of days due to difference in whole years days=0 if (Yr.gt.zeroYr) then do year=zeroYr,Yr-1 if (leapyr(year)) then days=days+366.0 else days=days+365.0 endif enddo endif * Now increment the number of days through the year days=days+real(MonthDays(Mo-1)+Dy-1) if (leapyr(Yr).and.(Mo.gt.2)) days=days+1.0 * Finally, the fraction of a day days=days+(real(Hr)+real(Mi)/60.0)/24.0 endif NumDays=days 101 format ((A),I) return end subroutine Date(zeroYr,days,Yr,Mo,Dy,Hr,Mi,MonthDays) * Given a number of days (since start of zeroYr), this routine works * out the true date implicit none integer zeroYr,Yr,Mo,Dy,Hr,Mi,MonthDays(0:12),dyear,Mday real*8 days,days1,Mdayprev logical fin,leapyr,ly * Calculate the number of years since start of zeroYr Yr=zeroYr days1=days fin=.false. 10 continue if (leapyr(Yr)) then dyear=366 else dyear=365 endif if(days1.ge.real(dyear)) then Yr=Yr+1 days1=days1-real(dyear) else fin=.true. endif if (.not.fin) goto 10 ly=leapyr(Yr) * Calculate the month number Mo=0 Mday=0 20 Mo=Mo+1 Mdayprev=Mday Mday=MonthDays(Mo) if (ly.and.(Mo.ge.2)) Mday=Mday+1 if (.not.(real(Mday).gt.days1)) goto 20 days1=days1-real(Mdayprev) * Calculate the day number Dy=int(days1)+1 days1=days1-real(Dy-1) * Calculate the hour print*,'days = ',days1 Hr=int(24.0*days1+0.01) * Calculate the minute Mi=INT(1440.0*days1-60.0*real(Hr)+0.5) return end logical function leapyr (yr) * Returns true do a leap year, false otherwise implicit none integer yr leapyr=((mod(yr,4).eq.0).and.((mod(yr,100).gt.0).or. : (mod(yr,400).eq.0))) return end subroutine PutInRange (angle,range,type) * To put angle in range -range to range if type=0 * 0 to 2*range if type=1 implicit none integer type real*8 angle,range,range2,mymod range2=2.*range if (type.eq.0) then angle = mymod(angle+range,range2)-range else angle = mymod(angle,range2) endif return end real*8 function mymod (x,interval) * This function puts x into the interval 0<=x 1 if (mainmin.eq.1) then * Minimum is leftmost one within bracket xxa=xa xxb=xa+real(indmins(1)+indmins(2))*(xb-xa)/200. * Store the fact that a minimum exists go=.true. elseif (mainmin.eq.nmins) then * Minimum is the rightmost one within bracket xxa=xa+real(indmins(nmins)+indmins(nmins-1))*(xb-xa)/200. xxb=xb * Store the fact that a minimum exists go=.true. else * The minimum is not the first or the last one xxa=xa+real(indmins(mainmin)+indmins(mainmin-1))*(xb-xa)/200. xxb=xa+real(indmins(mainmin)+indmins(mainmin+1))*(xb-xa)/200. * Store the fact that a minimum exists go=.true. endif else * Store the fact that no minima exist go=.false. endif if (go) then call GoldenMin1D(xxa,xxb,x,indx,y,y0,Error,Nobs,maxobs,aE,eE,iE, : Omega1E,omega2E,LE,Time,Epoch,long,lat,RefTime,tpe,obl) else if (verbose) then print*,'Note from subroutine PartialMin' print*,'No minima inside bracket ',xa,xb,' for parameter ', : indx endif endif if (verbose) then * Report partial J field and indicate position of minima print*,'----------------------------------------------' print*,'Partial min over index ',indx point=1 do l=0,100 if (l.eq.indmins(point)) then print 101,l,J(l),'*' if (point.lt.100) point=point+1 else print 102,l,J(l) endif enddo print*,'Main minimum is index ',indmins(mainmin) print*,'New value of parameter ',indx,' is ',x(indx) endif 101 format (I5,F12.4,(A)) 102 format (I5,F12.4) return end * ================================================================== subroutine GoldenMin1D (ax,cx,x,indx,y,y0,Error,Nobs,maxobs,aE,eE, : iE,Omega1E,omega2E,LE,Time,Epoch,long,lat, : RefTime,tpe,obl) * Golden ratio minimization in 1d * Requires function J(x) defined elsewhere * ------------------------------------------------------------------ * Parameters ax,cx (input: two x-values which bracket the min) * x(indx) (inout: location of minimum of J w.r.t. indx) * indx (input: the index of the parameter 1-6 varied) * y (passive: array to store predicted obs) * y0 (input: observations) * Error (input: observation errors) * Nobs (input: number of observations) * maxobs (input: maximum number of observations) * There must be a single minimum between ax and cx * ------------------------------------------------------------------ * Ross Bannister, April 2001 implicit none real*8 ratio,prec parameter (ratio=0.381966011,prec=0.0001) integer maxobs,loop,status,indx,Nobs,tpe(maxobs),l real*8 ax,cx,a,b,c,x(6),x0(6),CostFunc,delta,Jb,Jx,xx real*8 y(maxobs),y0(maxobs),Error(maxobs),RefTime,obl real*8 aE,eE,iE,Omega1E,omega2E,LE,Time(maxobs),Epoch,long,lat * Note arrangement of points (min between a and c) * a b c (status=0) * a b c (status=1) * Put the current x into the x0 structure do l=1,6 x0(l)=x(l) enddo a=ax c=cx loop=0 * Calculate intermediate point, b in first instance (put status=0) status=0 b=a+(c-a)*ratio 10 continue loop=loop+1 * Evaluate the cost function at point b x0(indx)=b call AstroFM (x0,aE,eE,iE,Omega1E,omega2E,LE,Time,Epoch, : long,lat,RefTime,tpe,Nobs,maxobs,y,obl) Jb=CostFunc(y,y0,Error,Nobs,maxobs) * Choose a new point, x and determine new bracketing and status if (status.eq.0) then xx=b+(c-b)*ratio * Evaluate the cost function at point xx x0(indx)=xx call AstroFM (x0,aE,eE,iE,Omega1E,omega2E,LE,Time,Epoch, : long,lat,RefTime,tpe,Nobs,maxobs,y,obl) Jx=CostFunc(y,y0,Error,Nobs,maxobs) if (Jb.lt.Jx) then * New arrangement is a,b,xx delta=c-xx c=xx status=1 else * New arrangement is b,xx,c delta=b-a a=b b=xx endif else xx=b-(b-a)*ratio x0(indx)=xx call AstroFM (x0,aE,eE,iE,Omega1E,omega2E,LE,Time,Epoch, : long,lat,RefTime,tpe,Nobs,maxobs,y,obl) Jx=CostFunc(y,y0,Error,Nobs,maxobs) if (Jb.lt.Jx) then * New arrangement is xx,b,c delta=xx-a a=xx status=0 else * New arrangement is a,xx,b delta=c-b c=b b=xx endif endif if (delta.gt.prec) goto 10 x(indx)=b return end * ================================================================== subroutine TimeRange (Time,Nobs,low,high,maxobs) * To return the lowest and highest time present in the time array implicit none integer Nobs,maxobs,loop,pairs,ind1 real*8 Time(maxobs),low,high pairs=Nobs/2 low=Time(1) high=low if (pairs.gt.1) then do loop=2,pairs ind1=(loop-1)*2+1 if (Time(ind1).gt.high) high=Time(ind1) if (Time(ind1).lt.low) low=Time(ind1) enddo endif return end * ================================================================== subroutine TrajCont (t1,t2,points,x,aE,eE,iE,Omega1E,omega2E,LE, : t0,obl,un) * To go over times between t1 and t2 (with 'points' intervals) and * find the RA/dec at each time (according to the model). implicit none integer points,loop,un real*8 t1,t2,x(6),aE,eE,iE,Omega1E,omega2E,LE,t0,t,ival,obl real*8 xcel,ycel,zcel,xecl,yecl,zecl,RA,dec,pi write (un,100) '#======================================='// : '==============================' write (un,100) '# '// : ' Ecliptic co-ordinates' write (un,100) '# Time (days) RA (model) Dec (model) '// : 'x (model) y (model) z(model)' write (un,100) '#---------------------------------------'// : '------------------------------' pi=4.*atan(1.) ival=(t2-t1)/real(points) do loop=1,points t=t1+real(loop-1)*ival * Find the gerocentric equatorial co-ords at this time call Orbit (x(1),x(2),x(3),x(4),x(5),x(6),aE,eE,iE,Omega1E, : omega2E,LE,t,t0,obl,xcel,ycel,zcel,xecl,yecl,zecl) * Convert to RA and Dec (radians) call Celestial (xcel,ycel,zcel,RA,dec) * Convert to hours,degrees RA=12.*RA/pi dec=180.*dec/pi write (un,101) t,RA,dec,xecl,yecl,zecl enddo write (un,100) '#======================================='// : '==============================' 100 format ((A)) 101 format (F15.6,2F12.6,3F10.3) return end * ================================================================== subroutine TrajPoint (Time,tpe,y0,Nobs,maxobs,reftime,long,lat,x, : aE,eE,iE,Omega1E,omega2E,LE,t0,obl,un) * Subroutine to output RA/dec obs (conversion is done if obs are of * type alt/azi) and projected obs at these times implicit none integer maxobs,tpe(maxobs),Nobs,loop,pairs,ind1,ind2,un real*8 Time(maxobs),y0(maxobs),long,lat,x(6),aE,eE,iE,Omega1E real*8 omega2E,LE,t0,obl,RA,dec,alt,azi,RAmodel,decmodel,altmodel real*8 azimodel,pi,reftime,deltaRA,deltaDec,deltaAlt,deltaAzi real*8 xcel,ycel,zcel,xecl,yecl,zecl,HalfCircle,HalfDay character*7 Label(2) Label(1)='RA/dec ' Label(2)='alt/azi' HalfCircle=180. HalfDay=12. write (un,100) : '#===========================================================' : //'===========================================================' write (un,100) : '# Time (days) RA(o) Dec(o) RA(m) Dec(m)' : //' RA(O-M)Dec(O-M) alt(o) azi(o) alt(m) ' : //' azi(m)alt(O-M)azi(O-M) Obs type' write (un,100) : '#-----------------------------------------------------------' : //'-----------------------------------------------------------' pi=4.*atan(1.0) pairs=Nobs/2 do loop=1,pairs print*,loop ind1=(loop-1)*2+1 ind2=ind1+1 * Deal with the observations if (tpe(ind1).eq.1) then * RA/dec variety of observation RA=12.*y0(ind1)/pi dec=180.*y0(ind2)/pi call local (cos(y0(ind2))*cos(y0(ind1)), : cos(y0(ind2))*sin(y0(ind1)), : sin(y0(ind2)),Time(ind1),long,lat,reftime,alt,azi) else * alt/azi variety of observation alt=180.*y0(ind1)/pi azi=180.*y0(ind2)/pi call localinverse (RA,dec,Time(ind1),long,lat,reftime,alt,azi) endif * Now make a prediction (model state) call Orbit (x(1),x(2),x(3),x(4),x(5),x(6),aE,eE,iE,Omega1E, : omega2E,LE,Time(ind1),t0,obl,xcel,ycel,zcel,xecl, : yecl,zecl) * Convert to RA and Dec (radians) call Celestial (xcel,ycel,zcel,RAmodel,decmodel) * Convert to hours,degrees RAmodel=12.*RAmodel/pi decmodel=180.*decmodel/pi * Convert to alt/azi call local (xcel,ycel,zcel,Time(ind1),long,lat,reftime,altmodel, : azimodel) altmodel=180.*altmodel/pi azimodel=180.*azimodel/pi * Compute residuals deltaRA=RA-RAmodel call PutInRange (deltaRA,HalfDay,0) deltaDec=dec-decmodel call PutInRange (deltaDec,HalfCircle,0) deltaAlt=alt-altmodel call PutInRange (deltaAlt,HalfCircle,0) deltaAzi=azi-azimodel call PutInRange (deltaAzi,HalfCircle,0) * Now output the data write (un,101) Time(ind1),RA,dec,RAmodel,decmodel, : deltaRA,deltaDec,alt,azi,altmodel, : azimodel,deltaAlt,DeltaAzi, : ' ',Label(tpe(ind1)) enddo write (un,100) : '#===========================================================' : //'===========================================================' 100 format ((A)) 101 format (F15.6,12F8.2,A2,A7) return end * ================================================================== subroutine SecondH (x0,aE,eE,iE,Omega1E,omega2E,LE,Time,t0, : long,lat,reftime,tpe,Nobs,maxobs,Hd2,obl, : ObsA,ObsB,ObsC,ObsD) * Compute the maxobsx6x6 matrix of second derivatives of the * forward model * December 2002, Ross Bannister implicit none integer maxobs,tpe(maxobs),Nobs,i,j,ob real*8 x0(6),aE,eE,iE,Omega1E,omega2E,LE,Time(maxobs),t0,long,lat real*8 reftime,Hd2(maxobs,6,6),obl real*8 delta,xA(6),xB(6),xC(6),xD(6) real*8 ObsA(maxobs),ObsB(maxobs),ObsC(maxobs),ObsD(maxobs) do i=1,6 xA(i)=x0(i) xB(i)=x0(i) xC(i)=x0(i) xD(i)=x0(i) enddo delta=0.05 * Compute the second derivatives of the forward model * Exploit symmetry of second derivatives do i=1,6 do j=1,6 if (j.gt.i) then * Upper-triangular components * Set-up the states for computation of the second derivatives xA(i)=x0(i)+delta xA(j)=x0(j)+delta xB(i)=x0(i)-delta xB(j)=x0(j)+delta xC(i)=x0(i)+delta xC(j)=x0(j)-delta xD(i)=x0(i)-delta xD(j)=x0(j)-delta * Call the forward model with these incremented state vectors call AstroFM (xA,aE,eE,iE,Omega1E,omega2E,LE,Time,t0, : long,lat,reftime,tpe,Nobs,maxobs,ObsA,obl) call AstroFM (xB,aE,eE,iE,Omega1E,omega2E,LE,Time,t0, : long,lat,reftime,tpe,Nobs,maxobs,ObsB,obl) call AstroFM (xC,aE,eE,iE,Omega1E,omega2E,LE,Time,t0, : long,lat,reftime,tpe,Nobs,maxobs,ObsC,obl) call AstroFM (xD,aE,eE,iE,Omega1E,omega2E,LE,Time,t0, : long,lat,reftime,tpe,Nobs,maxobs,ObsD,obl) * Compute the second derivatives do ob=1,Nobs Hd2(ob,i,j)=(ObsA(ob)-ObsB(ob)-ObsC(ob)+ObsD(ob))/ : (4.*delta*delta) enddo * Return states to original values xA(i)=x0(i) xA(j)=x0(j) xB(i)=x0(i) xB(j)=x0(j) xC(i)=x0(i) xC(j)=x0(j) xD(i)=x0(i) xD(j)=x0(j) elseif (j.eq.i) then * Diagonal components * Set-up the states for computation of the second derivatives xA(i)=x0(i)+delta xB(i)=x0(i)-delta * Call the forward model with these incremented state vectors call AstroFM (xA,aE,eE,iE,Omega1E,omega2E,LE,Time,t0, : long,lat,reftime,tpe,Nobs,maxobs,ObsA,obl) call AstroFM (xB,aE,eE,iE,Omega1E,omega2E,LE,Time,t0, : long,lat,reftime,tpe,Nobs,maxobs,ObsB,obl) call AstroFM (x0,aE,eE,iE,Omega1E,omega2E,LE,Time,t0, : long,lat,reftime,tpe,Nobs,maxobs,ObsC,obl) * Compute the second derivatives do ob=1,Nobs Hd2(ob,i,j)=(ObsA(ob)-2.*ObsC(ob)+ObsB(ob))/(delta*delta) enddo * Return states to original values xA(i)=x0(i) xB(i)=x0(i) else * Exploit symmetry of second derivatives do ob=1,Nobs Hd2(ob,i,j)=Hd2(ob,j,i) enddo endif enddo enddo return end * ================================================================== subroutine GaussEll (A,A1,b,b1,x,N,Nmax) * Solve for x in Ax = b (matrix A vector b) (Gaussian ellimination) * Ross Bannister, December 2002 implicit none integer N,Nmax,i,ip,j real*8 A(Nmax,Nmax),A1(Nmax,Nmax),b(Nmax),b1(Nmax),x(Nmax),factor real*8 sum * Copy A into A1 and b into b1 (workspace arrays) do i=1,N do j=1,N A1(i,j)=A(i,j) enddo b1(i)=b(i) enddo * Stage 1 of the solver : preparation to upper triangular form do i=1,N-1 do ip=i+1,N factor=A1(ip,i)/A1(i,i) * Do replacements do j=i+1,N A1(ip,j)=A1(ip,j)-factor*A1(i,j) enddo b1(ip)=b1(ip)-factor*b1(i) * Elliminate terms A1(ip,i)=0. enddo enddo * Stage 2 of the solver : find x x(N)=b1(N)/A1(N,N) do i=N-1,1,-1 sum=0. do j=i+1,N sum=sum+A1(i,j)*x(j) enddo x(i)=(b1(i)-sum)/A1(i,i) enddo * The correct answer should now be in x return end * ================================================================== subroutine Jacobi (A,B,T,N,max,lim,trace) * To diagonalize a square symmetric matrix by a sequence of Jacobi * rotations * Ross Bannister, January 2001 / August 2002 * A (input) matrix in original representation * B (output) matrix in diagonal representation * T (output) matrix of eigenvectors (rows) * N (input) the order of the matrix * max (input) the maximum order of the matrix (for array dim) * lim (input) convergence criteria, e.g. lim=0.01 * trace (output) Sum of eigenvalues implicit none integer N,max,i,j,p,q,k,sweeps * integer lp1,lp2 real*8 A(max,max),B(max,max),T(max,max),lim,off,diag,theta real*8 t1,t2,tt,cc,ss,sc,Bip,Bqj,Tpj,Bpp,Bqq,dis,s,c,norm,trace if (N.gt.max) then print*,'Report from matrix JACOBI.' print*,'N =',N print*,'max =',max print*,'Please increase the value of max.' stop endif * The matrix B is initially equal to A do j=1,N do i=1,N B(i,j)=A(i,j) enddo enddo * The matrix T is initially the identity matrix do j=1,N do i=1,N if (i.eq.j) then T(i,j)=1.0 else T(i,j)=0.0 endif enddo enddo * Calculate the sum of the square of off-diagonal elements of B off=0.0 do j=2,N do i=1,j-1 off=off+B(i,j)*B(i,j) enddo enddo off=2.0*off * Calculate the sum of the square of diagonal elements of B diag=0.0 do i=1,N diag=diag+B(i,i)*B(i,i) enddo * The number of sweeps completed sweeps=0 10 continue * Loop round each upper-right element (a 'sweep') do p=2,N do q=1,p-1 if (abs(B(q,p)).gt.(0.0000001)) then * In order to find the new representation of the matrix, * work out some details regarding the rotation, U to eliminate * element q,p of B (p>q). theta=(B(p,p)-B(q,q))/(2.0*B(q,p)) dis=sqrt(theta*theta+1) * Two possible values of t t1=-theta+dis t2=-theta-dis * Choose the smaller value if (abs(t1).lt.abs(t2)) then tt=t1 else tt=t2 endif * Find the values of s and c c=1/sqrt(tt*tt+1) s=tt*c * Modify the operator UBU(trans) (piecewise) * Update only upper-right part of matrix cc=c*c ss=s*s sc=s*c * Columns p and q: do i=1,p-1 if (i.ne.q) then Bip=B(i,p) B(i,p)=c*Bip+s*B(i,q) if (i.le.q) B(i,q)=c*B(i,q)-s*Bip endif enddo * Rows p and q: do j=q+1,N if (j.ne.p) then Bqj=B(q,j) B(q,j)=c*Bqj-s*B(p,j) if (j.gt.p) B(p,j)=c*B(p,j)+s*Bqj endif enddo * Mixed elements * B(q,p)=(cc-ss)*B(q,p)+sc*(B(q,q)-B(p,p)) B(q,p)=0.0 * New diagonal elements Bpp=B(p,p) Bqq=B(q,q) B(p,p)=cc*Bpp+ss*B(q,q)+2.0*sc*B(p,q) B(q,q)=cc*B(q,q)+ss*Bpp-2.0*sc*B(p,q) * Modified value of sum of square of off-diagonal elements... off=off-2.0*B(p,q)*B(p,q) * ...and diagonal elements diag=diag-Bpp*Bpp-Bqq*Bqq+B(p,p)*B(p,p)+B(q,q)*B(q,q) * Update lower-left of matrix do k=1,N if (k.lt.p) then B(p,k)=B(k,p) if (k.lt.q) B(q,k)=B(k,q) endif if (k.gt.q) then B(k,q)=B(q,k) if (k.gt.p) B(k,p)=B(p,k) endif enddo * Modify the transformation T(new)=U T(old) * Only rows p and q are affected do j=1,N Tpj=T(p,j) T(p,j)=c*Tpj+s*T(q,j) T(q,j)=c*T(q,j)-s*Tpj enddo endif enddo enddo sweeps=sweeps+1 if (((off/diag).gt.lim).or.(sweeps.lt.3)) goto 10 * Normalize the eigenvectors do i=1,N norm=0.0 do j=1,N norm=norm+T(i,j)*T(i,j) enddo norm=sqrt(norm) do j=1,N T(i,j)=T(i,j)/norm enddo enddo * Calculate the sum of eigenvalues trace=0. do i=1,6 trace=trace+B(i,i) enddo return end * ================================================================== real*8 function son (input) * Returns the sign of the input +1.,-1.,0. implicit none real*8 input,result if (input.gt.(0.)) then result=1. elseif (input.lt.(0.)) then result=-1. else result=0. endif son=result return end * ==================================================================