c SUBROUTINE TO COMPUTE THE TRAVEL TIMES c BASED ON THE DEFINITION c BY GIZON & BIRCH (2004) EQUATIONS 3 & 4 c VERSION 1.0, 03/29/2007 c AUTHOR SEBASTIEN COUVIDAT c ISSUES TO SOLVE c the format of the cross-covariance file has not been discussed c the width and center of time windows has not been discussed c the name and format of the output file have not been discussed program gbtimes implicit none double precision pi,dt integer nx,ny,nt,naxis,nmax parameter (pi=3.141592653589793116d0) parameter (naxis=3,nmax=4096) parameter (nx=256,ny=256,nt=127,dt=60.d0) integer status,unit,readwrite,naxes(naxis),iii,i,j,k,blocksize,nfound double precision correl(nx,ny,nt),Weightp(nt),Weightm(nt), 1 position(nt),denominatorp,denominatorm,windowp(nt),taup(nx,ny), 2 taum(nx,ny),windowm(nt),correlref(nt),correlrefder(nt),tau(nx,ny,2) character*100 filename,filewrite external readfits,writefits,ftgiou,ftopen,ftgknj c open the FITS file containing the cross-covariance filename='correlation_MidRes_im3dV2_shifted_5.fits' !format correl[256,256,127] negative, then positive times filewrite='test.fits' iii=27 status=0 call ftgiou(unit,status) readwrite=0 call ftopen(unit,filename,readwrite,blocksize,status) call ftgknj(unit,'NAXIS',1,naxis,naxes,nfound,status) if (nfound .ne. 3)then print *,'READIMAGE failed to read the NAXISn keywords.' endif if ((naxes(1) .ne. nx) .OR. (naxes(2) .ne. ny) .OR. (naxes(3) .ne. nt))then print *,'PROBLEM: DIMENSIONS NOT CORRECT' endif call readfits(unit,nx,ny,nt,correl) call window(windowp,windowm,nt,iii,dt) c compute the reference cross-covariance do k=1,nt correlref(k) = 0.d0 do i=1,nx do j=1,ny correlref(k)=correlref(k)+correl(i,j,k) enddo enddo correlref(k)=correlref(k)/dble(nx)/dble(ny) enddo do k=1,nt print *,(k-1.d0)*dt,windowp(k),windowm(k) enddo c compute the time derivative of the reference cross-covariance c call differentation(correlref,dt,nt,correlrefder) call differentation2(correlref,nt,dt,correlrefder) c compute the weight functions denominatorp = 0.d0 do k=1,nt denominatorp = denominatorp+dt*windowp(k)*correlrefder(k)**2 enddo denominatorm = 0.d0 do k=1,nt denominatorm = denominatorm+dt*windowm(k)*correlrefder(k)**2 enddo do k=1,nt Weightp(k)= -1.d0*windowp(k)*correlrefder(k)/denominatorp Weightm(k)= windowm(k)*correlrefder(k)/denominatorm enddo c compute the travel times do i=1,nx do j=1,nx taup(i,j) = 0.d0 taum(i,j) = 0.d0 do k=1,nt taup(i,j) = taup(i,j)+dt*Weightp(k)*(correl(i,j,k)-correlref(k)) taum(i,j) = taum(i,j)+dt*Weightm(k)*(correl(i,j,k)-correlref(k)) enddo tau(i,j,1) = (taup(i,j)+taum(i,j))/2.d0 tau(i,j,2) = taup(i,j)-taum(i,j) enddo enddo c write the result file call writefits(filewrite,nx,ny,naxis,tau) stop end c ------------------------------------------------------------------------------- c computes a derivative subroutine differentation(f,h,n,fp) double precision f(n),h,fp(n) fp(1) = 1.d0/(12.d0*h)*(-25.d0*f(1) + 48.d0*f(2) - 36.d0*f(3) + 16.d0*f(4) - 3.d0*f(5)) fp(2) = 1.d0/(12.d0*h)*(-3.d0*f(1) - 10.d0*f(2) + 18.d0*f(3) - 6.d0*f(4) + f(5)) do i=3,n-2 fp(i) = 1.d0/(12.d0*h)*(f(i-2) - 8.d0*f(i-1) + 8.d0*f(i+1) - f(i+2) ) enddo fp(nt-1)= 1.d0/(12.d0*h)*(-f(nt-4) + 6.d0*f(nt-3) - 18.d0*f(nt-2) + 10.d0*f(nt-1) + 3.d0*f(nt) ) fp(nt) = 1.d0/(12.d0*h)*(3.d0*f(nt-4) - 16.d0*f(nt-3) + 36.d0*f(nt-2) - 48.d0*f(nt-1) + 25.d0*f(nt)) return end c does a Fourier interpolaiton by a factor 16, and then computes the derivative subroutine differentation2(f,nt0,dt0,fp0) parameter (pi=3.141592653589793116d0,factor=64.d0,nt=8128) ! nt=factor*nt0 integer nt0 double precision f(nt0),fi(nt),var,varr,time(nt),time0(nt0) double precision fp(nt),h,fp0(nt0),dt0,dt dt = dt0/factor do i=1,nt time(i)=(i-1.d0)*dt enddo do i=1,nt0 time0(i)=(i-1.d0)*dt0 enddo do i=1,nt var=0.d0 do j=1,nt0 if ((time(i)-time0(j)) .eq. 0.d0) then varr=f(j) else varr=f(j)*sin(pi*(time(i)-time0(j))/dt0)/(time(i)-time0(j))*dt0/pi endif var=var+varr enddo fi(i) = var enddo fp(1) = 1.d0/(12.d0*dt)*(-25.d0*fi(1) + 48.d0*fi(2) - 36.d0*fi(3) + 16.d0*fi(4) - 3.d0*fi(5)) fp(2) = 1.d0/(12.d0*dt)*(-3.d0*fi(1) - 10.d0*fi(2) + 18.d0*fi(3) - 6.d0*fi(4) + fi(5)) do i=3,nt-2 fp(i) = 1.d0/(12.d0*dt)*(fi(i-2) - 8.d0*fi(i-1) + 8.d0*fi(i+1) - fi(i+2) ) enddo fp(nt-1)= 1.d0/(12.d0*dt)*(-fi(nt-4) + 6.d0*fi(nt-3) - 18.d0*fi(nt-2) + 10.d0*fi(nt-1) + 3.d0*fi(nt) ) fp(nt) = 1.d0/(12.d0*dt)*(3.d0*fi(nt-4) - 16.d0*fi(nt-3) + 36.d0*fi(nt-2) - 48.d0*fi(nt-1) + 25.d0*fi(nt)) do i=1,nt0 fp0(i)=fp((i-1)*nt/nt0+1) enddo return end c ------------------------------------------------------------------------------- subroutine window(windowp,windowm,nt,iii,dt) double precision windowp(nt),windowm(nt),centrage,dt,time(nt) integer iii,k,nt do k=1,nt time(k) = (k-1.d0)*dt-(nt+1.d0)/2.d0*dt ! for nt odd enddo if (iii .EQ. 2) then centrage=1140.d0 else if (iii .EQ. 7) then centrage=1150.d0 else if (iii .EQ. 12) then centrage=1200.d0 else if (iii .EQ. 17) then centrage=1500.d0 else if (iii .EQ. 22) then centrage=1650.d0 else if (iii .EQ. 27) then centrage=1750.d0 else if (iii .EQ. 32) then centrage=1850.d0 else if (iii .EQ. 37) then centrage=2000.d0 else if (iii .EQ. 42) then centrage=2100.d0 else if (iii .EQ. 47) then centrage=2200.d0 else if (iii .EQ. 52) then centrage=2300.d0 endif print *,'CENTER OF TIME WINDOW',centrage do k=1,nt if ((ABS(time(k)-centrage) .le. 600) .AND. (k .ge.(nt+1.d0)/2.d0)) then windowp(k)=1.d0 else windowp(k)=0.d0 endif if ((ABS(time(k)+centrage) .le. 600) .AND. (k .lt.(nt+1.d0)/2.d0)) then windowm(k)=1.d0 else windowm(k)=0.d0 endif enddo return end c -------------------------------------------------------------------------------- subroutine printerror(status) C This subroutine prints out the descriptive text corresponding to the C error status value and prints out the contents of the internal C error message stack generated by FITSIO whenever an error occurs. integer status character errtext*30,errmessage*80 C Check if status is OK (no error); if so, simply return if (status .le. 0)return C The FTGERR subroutine returns a descriptive 30-character text string that C corresponds to the integer error status number. A complete list of all C the error numbers can be found in the back of the FITSIO User's Guide. call ftgerr(status,errtext) print *,'FITSIO Error Status =',status,': ',errtext C FITSIO usually generates an internal stack of error messages whenever C an error occurs. These messages provide much more information on the C cause of the problem than can be provided by the single integer error C status value. The FTGMSG subroutine retrieves the oldest message from C the stack and shifts any remaining messages on the stack down one C position. FTGMSG is called repeatedly until a blank message is C returned, which indicates that the stack is empty. Each error message C may be up to 80 characters in length. Another subroutine, called C FTCMSG, is available to simply clear the whole error message stack in C cases where one is not interested in the contents. call ftgmsg(errmessage) do while (errmessage .ne. ' ') print *,errmessage call ftgmsg(errmessage) end do end c --------------------------------------------------------------------------- subroutine deletefile(filename,status) C A simple little routine to delete a FITS file integer status,unit,blocksize character*(*) filename C Simply return if status is greater than zero if (status .gt. 0)return C Get an unused Logical Unit Number to use to open the FITS file call ftgiou(unit,status) C Try to open the file, to see if it exists call ftopen(unit,filename,1,blocksize,status) if (status .eq. 0)then C file was opened; so now delete it call ftdelt(unit,status) else if (status .eq. 103)then C file doesn't exist, so just reset status to zero and clear errors status=0 call ftcmsg else C there was some other error opening the file; delete the file anyway status=0 call ftcmsg call ftdelt(unit,status) end if C Free the unit number for later reuse call ftfiou(unit, status) end c ----------------------------------------------------------------------------------- subroutine writefits(filename,n1,n2,naxis,array) integer status,unit,blocksize,bitpix,naxis,naxes(naxis) integer i,j,group,fpixel,nelements,n1,n2 character filename*80 logical simple,extend double precision array(0:(n1-1),0:(n2-1),2) status=0 call deletefile(filename,status) call ftgiou(unit,status) blocksize=1 call ftinit(unit,filename,blocksize,status) simple=.true. bitpix=-32 c naxis=2 naxes(1)=n1 naxes(2)=n2 naxes(3)=2 nelements=n1*n2*2 extend=.true. call ftphpr(unit,simple,bitpix,naxis,naxes,0,1,extend,status) group=1 fpixel=1 c nelements=naxes(1)*naxes(2)*naxes(3) c nelements=naxes(1)*naxes(2) call ftpprd(unit,group,fpixel,nelements,array,status) call ftclos(unit, status) call ftfiou(unit, status) if (status .gt. 0)call printerror(status) end c --------------------------------------------------------------------------------- subroutine readfits(unit,nx,ny,nt,array) integer status,unit,nx,ny,nt,ix,iy,it integer group,firstpix,nbuffer,npixels real nullval double precision array(0:nx-1,0:ny-1,0:nt-1) c real*8 array(0:nx-1,0:ny-1,0:nt-1) double precision buffer(0:nx-1) c real*8 buffer(0:nx-1),array(0:(2*nx)-1,0:(2*ny)-1,0:(2*nt)-1) c real*8 buffer(0:nx-1),array(0:nx-1,0:ny-1,0:nt-1) c real allocatable::buffer(:),array(:,:,:) logical anynull character filename*80 status=0 C Initialize variables npixels=nx*ny*nt group=1 firstpix=1 nullval=0 nbuffer=nx do it=0,nt-1 do iy=0,ny-1 call ftgpvd(unit,group,firstpix,nbuffer,nullval, 1 buffer,anynull,status) firstpix=firstpix + nbuffer do ix=0,nx-1 array(ix,iy,it)=buffer(ix) enddo enddo enddo C The FITS file must always be closed before exiting the program. C Any unit numbers allocated with FTGIOU must be freed with FTFIOU. call ftclos(unit, status) call ftfiou(unit, status) C Check for any error, and if so print out error messages. C The PRINTERROR subroutine is listed near the end of this file. if (status .gt. 0)call printerror(status) return end c ------------------------------------------------------------------------