c     This program corrects the shot table for the gun delay
c     and for the offset between GPS and airguns.
c     In addition, shotpoint numbers are added.
c
c     CCGS Hudson  June 24, 2001
c
      character*60 input,output
      REAL K0,FE,NORTH
	real*8 slat,slong
      INTEGER*2 calday1,calhr1,calmin1,calday2,calhr2,calmin2,
     *          gundelay,clockoff
	INTEGER*2 startday,starthr,startmin,endday,endhr,endmin

      COMMON /CONST/A0,A2,A4,A6,A8
      COMMON /CONST2/A,E1,E2
C
C     ELLIPSIODAL PARAMETERS FROM   CLARKE 1866 ELLIPSIOD
C 
      A=6378206.4
      RECIPF=294.97898
      K0=0.9996
      RADDEG=57.29577951
      FE=500000.
      F=1./RECIPF
      B=A-A*F
      E1=(A*A-B*B)/(A*A)
      E2=(A*A-B*B)/(B*B)
C     ARC LENGTH CONSTANTS
      A0=1.-0.25*E1-0.046875*E1*E1-0.01953125*E1**3.-175./16384.*E1**4.
      A2=0.375*(E1+0.25*E1*E1+0.1171875*E1**3.-455./4096.*E1**4.)
      A4=0.05859375*(E1*E1+0.75*E1**3.-77./128.*E1**4.)
      A6=35./3072.*(E1**3.-1.28125*E1**4.)
      A8=-315./131072.*E1**4.
C
c
      pi=3.141592653589
c
c-------------------------------------------------------------
c
c     Distance between GPS antenna and midship line in meter
c     (antenna on port side = minus, starboard = +)
c
      ds=3.3
c
c     Distance of airguns behind GPS antenna (in meter)
c
      da=58.1+30.5
c
c     Delay of airgun shots (in sec)
c
      dt=0.000
c
c--------------------------------------------------------------
c
c     Opening IO files
c
      write (*,*) 'Number of first shot?'
      read (*,*) ish
c
      write (*,*) 'Line number?'
      read (*,*) iline
c
   15 write (*,*) 'Name of input file:'
      read (*,999) input
  999 format (a60) 
      open (10,file=input,status='old',err=10)
      goto 20
   10 write (*,*) 'File does not exist. Try again.'
      goto 15
c
   20 write (*,*) 'Name of output file:'
      read (*,999) output
      open (11,file=output,status='new',err=30)
      goto 40
   30 write (*,*) 'File exists already. Try again.'
      goto 20
c
c--------------------------------------------------------------
c
c     Reading raw navigation
c
 40   write (11,996) iline
 996  format ('Scotian Margin 2001 - Mariprobe Line ',i1)
c
 45   read (10,998,err=50,end=100) id,ih,im,sec,x,y,h
 998  format (i4,i3,1x,i2,1x,f6.3,2f11.5,f6.1)
c
c     Set central meridian
c
      cm=x
c
c--------------------------------------------------------------
c
c     Apply time correction
c
      sec=sec+dt
      if (sec.ge.60.) then
          sec=sec-60.
          im=im+1
          if (im.ge.60) then
              im=im-60
              ih=ih+1
              if (ih.ge.24) then
                  ih=ih-24
                  id=id+1
                else
              endif
            else
          endif
        else
      endif
c
c--------------------------------------------------------------
c
c     Conversion of geographical coordinates to UTM
c
      call llutm (y,x,cm,uy,ux)
c
c--------------------------------------------------------------
c
c     Correction for lateral offset of GPS antenna from
c     midship line
c
      dsy=ds*sin(pi*h/180.)
      dsx=-ds*cos(pi*h/180.)
      ux=ux+dsx
      uy=uy+dsy
c
c--------------------------------------------------------------
c
c     Correction for airgun offset (behind antenna)
c
      day=-da*cos(pi*h/180.)
      dax=-da*sin(pi*h/180.)
      ux=ux+dax
      uy=uy+day
c
c--------------------------------------------------------------
c
c     Back projection from UTM to geographical coordinates
c
      call utmll (y,x,uy,ux,cm)
c
c--------------------------------------------------------------
c
c     Writing output file
c
      write (11,997) ish,id,ih,im,sec,x,y,h
 997  format (i5,i4,i3,':',i2,':',f6.3,2f11.5,f6.1)
      ish=ish+1
      goto 45
c
c--------------------------------------------------------------
c
c     Ending program
c
 50   write (*,*) 'Error in reading input file.'
 100  close (10)
      close (11)
      stop
      end
c
c
c
C**************************************************
C
      SUBROUTINE UTMLL(LAT,LON,NORTH,EAST,CM)
 
C     COMPUTE LATITUDE/LONGITUDE FROM UTMS

      REAL RADDEG,DEGRAD,NORTH,EAST,CM,PHIFT,T,M,N,N1,T1,T2,T3,LAT,LON
     +     ,FE,K0
      COMMON /CONST2/ A,E1,E2
      DATA DEGRAD/1.7453292519943E-2/,FE/500000./,K0/0.9996/
 
      RADDEG=1./DEGRAD
      CALL FOOTPT(NORTH,PHIFT)
      PHIFT=PHIFT*DEGRAD
      Y=NORTH/K0
      X=(EAST-FE)/K0
      T=TAN(PHIFT)
C  MERIDIAN RADIUS OF CURVATURE AT FOOT POINT LATITUDE
      M=A*(1-E1)/(1-E1*SIN(PHIFT)**2.)**1.5
C  PRIME VERTICAL RADIUS OF CURVATURE AT FOOT POINT LATITUDE
      N1=A/SQRT(1-E1*SIN(PHIFT)**2.)
      N=E2*COS(PHIFT)**2.
      T1=-T*X/M*X/N1/2. + T*X/M*X/N1*X/N1*X/N1/24.*(5.+3.*T*T+N+4.*N*N+
     +9.*N*T*T)
      T2=T*X/M*X/N1*X/N1*X/N1*X/N1*X/N1/720.*(61.-90.*T*T+46.*N+45.*
     +T**4.-252.*T*T*N+3.*N*N+100.*N**3.-66.*T*T*N*N-90.*T**4.*N+88.*
     +N**4.+225.*T**4.*N*N+84.*T*T*N**3.-192.*T*T*N**4.)
      T3=T*X/M*X/N1*X/N1*X/N1*X/N1*X/N1*X/N1*X/N1/40320.*(1385.+3633.*
     +T*T+4095.*T**4+1575.*T**6)
      LAT=(PHIFT+T1-T2+T3)*RADDEG
C
      T1=X/N1-X/N1*X/N1*X/N1/6.*(1.+2.*T*T+N)
      T2=5.+6.*N+28.*T*T-3.*N*N+8.*T*T*N+24.*T**4.-4.*N**3.+4.*T*T*N*N+
     +24.*T*T*N**3
      T3=61.+662.*T*T+1320.*T**4.+720.*T**6.
      LON=CM+(T1+X/N1*X/N1*X/N1*X/N1*X/N1/120.*T2-X/N1*X/N1*X/N1*X/N1*
     +X/N1*X/N1*X/N1/5040.*T3)/COS(PHIFT)*RADDEG
      RETURN
      END
 
C
C******************************************
C
      SUBROUTINE FOOTPT(NORTH,PHIFT)
 
      REAL LATFT,PHI,NORTH,PHIFT,MEDARC,K0
 
C   CALCULATES FOOT POINT LATITUDE

      DATA DEGRAD/1.7453292519945E-2/,K0/0.9996/

      LATFT=INT(NORTH*90./10000000.)+1.
      PHI=(LATFT-1.)*DEGRAD
      ARC=MEDARC(PHI)*K0
 200  IF (NORTH.GE.ARC) GOTO 100
        LATFT=LATFT-1.
        PHI=LATFT*DEGRAD
        ARC=MEDARC(PHI)*K0
        GOTO 200
 100  LATFT=LATFT-1.
      PHIFT=LATFT
      PHI=LATFT*DEGRAD
      ARC=MEDARC(PHI)*K0
      LATFT=LATFT+1.
      PHI=LATFT*DEGRAD
      ARC2=MEDARC(PHI)*K0
      PHIFT=PHIFT+(NORTH-ARC)/(ARC2-ARC)
      PHI=PHIFT*DEGRAD
      ARC3=MEDARC(PHI)*K0
      PHIFT=PHIFT+(NORTH-ARC3)/(ARC2-ARC)
      RETURN
      END
C
C**************************************
C
 
      REAL FUNCTION MEDARC(PHI)
 
      REAL  PHI
      COMMON /CONST/A0,A2,A4,A6,A8
      COMMON /CONST2/A,E1,E2
 
      B=A0*PHI-A2*SIN(2*PHI)+A4*SIN(4*PHI)-A6*SIN(6*PHI)+A8*SIN(8*PHI)
      MEDARC=A*B
      RETURN
      END
 
C
C*********************************
C
      SUBROUTINE LLUTM(LAT,LON,CM,NORTH,EAST)
C
C  THIS SUBROUTINE CONVERTS DECIMAT LATITUDES AND LONGITUDES TO
C  STANDARD U.T.M. COORDINATES
C
      REAL LAT,LON,NORTH,EAST,N,N2,FE,K0,MEDARC
      COMMON /CONST2/A,E1,E2
 
      PHI=LAT/57.29577951
      V=E1*SIN(PHI)**2.
      N=A/SQRT((1.-V))
      T=TAN(PHI)**2.
      DELTA=(LON-CM)/57.29577951
      N2=E2*COS(PHI)**2.
      K0=0.9996
      S0=MEDARC(PHI)
      FE=0.
 
      T1=DELTA*DELTA/2.*SIN(PHI)*COS(PHI)
      T2=DELTA**4./24.* SIN(PHI)*COS(PHI)**3.*(5.-T+9.*N2+4.*N2*N2)
      T3=61.-58.*T+T*T+270.*N2-330.*T*N2+445.*N2*N2+324.*N2**3.-680.*
     + N2*N2*T+88.*N2**4.-600.*N2**3.*T-192.*N2**4.*T
      T4=1385.-3111.*T+543.*T*T-T**3.
      NORTH=K0*S0+K0*N*(T1+T2+DELTA**6./720.*SIN(PHI)*COS(PHI)**5.*T3+
     + DELTA**8./40320.*SIN(PHI)*COS(PHI)**7.*T4)
      T1=DELTA*COS(PHI)+(DELTA*COS(PHI))**3./6.*(1-T+N2)
      T2=5.-18.*T+T*T+14.*N2-58.*T*N2+13.*N2*N2+4.*N2**3.-64.*N2*N2*T-
     + 24.*N2**3.*T
      T3=(61.-479.*T+179*T*T-T**3.)/5040.
      EAST=K0*N*(T1+(DELTA*COS(PHI))**5./120.*T2+(DELTA*COS(PHI))**7.
     + *T3)
      EAST=500000.+EAST
      RETURN
      END
