      PROGRAM SHOTUTM
C
C     This version reads shot table file and replace its inacurate
c     position with the position of ships nav data (incl. ship's heading
c     either raw or calculated using head.for)
c     Revised: 26 Nov 1995 to include correction for offset between
c     GPS position and location of guns (shtcor in metres)
c	This version for input of SHIPS Navigation File where lat,long
c	refers to position of guns.

C      UTM <--> LAT/LON ROUTINES
C      FROM FORMULAE PROVIDED BY UNB SCHOOL OF SURVEYING
C      FROM SURVEYS AND MAPPING DIVISION OF THE
C      LAND REGISTRATION AND INFORMATION SERVICE
C

      REAL K0,FE,NORTH
	real*8 slat,slong
      CHARACTER ID*80,KIN1*40,KOUT*40
      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
      DATA shtcor/0./
C
C     ELLIPSIODAL PARAMETERS FROM   CLARKE 1866 ELLIPSIOD
 
      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     INITIALISATION
C
      PRINT 100
  100 FORMAT(///,' PROGRAM SHOTFORM  ',///)
  15  PRINT 101
  101 FORMAT(/,' WHAT IS THE NAME OF THE GPSTAR INPUT SHOTLOG FILE? ')
      READ(*,102) KIN1
  102 FORMAT(A40)
C
      OPEN(2,FILE=KIN1,STATUS='OLD',FORM='FORMATTED')
C
      print *, 'Enter offset in metres between GPS antenna and gun '
	read(*,*) shtcor
c
  20  PRINT 103
  103 FORMAT(/,' WHAT IS THE NAME OF THE OUTPUT SHOTFILE FOR DREA2SGY?')
      READ(*,102) KOUT
      OPEN(3,FILE=KOUT,ERR=21,STATUS='UNKNOWN',FORM='FORMATTED')
      GOTO 22
  21  PRINT *, ' ERROR OPENING FILE ',KOUT
      GOTO 20
  22  CONTINUE
C
      READ(2,'(A80)') ID
C
      PRINT *, ' ENTER DECIMAL LAT, LONG OF OBS, AND CENTRAL MERIDIAN: '
      READ *, OBSLAT,OBSLON,CM
      print *, ' Enter time of clock reset (dayno,hr,min,sec): '
      read *, calday1,calhr1,calmin1,calsec1
      print *, ' Enter time of drift calibration (dayno,hr,min,sec): '
      read *, calday2,calhr2,calmin2,calsec2
      print *, 'Enter total clock drift (in msec)'
      print *, '   +/- mean that obs is fast/slow: '
      read *, drift
      print *,' Enter gun delay (in msec): '
      read *, gundelay
      print *,' Enter offset between ship and obs clocks (in msec)'
      print *,'   +/- mean that obs is fast/slow w/r/t ship: '
      read *, clockoff
c
	read *, startday,starthr,startmin
	read *, endday,endhr,endmin
	stime=startday*24+starthr+startmin/60.
	etime=endday*24+endhr+endmin/60.
c
c      stime=1.
c      etime=8000.
c
      t1=calday1*24.+calhr1+calmin1/60.+calsec1/3600.
      t2=calday2*24.+calhr2+calmin2/60.+calsec2/3600.
      delt=t2-t1
      driftrate=drift/delt

C      CLATOBS=COS(obslat/RADDEG)
C      SLATOBS=SIN(obslat/RADDEG)

      CALL LLUTM(OBSLAT,OBSLON,CM,OBSNORTH,OBSEAST)
C      print *, obsnorth,obseast
C      pause

C     print header

      write(3,'(a80)') id
      write(3,"(2f12.6,2X,'OBS Ref. latitude and longitude')")
     *         obslat,obslon
      write(3,"(3i3,f4.0,9x,'Time of clock calibration')")
     *         calday1,calhr1,calmin1,calsec1
      write(3,"(f10.4,12x,'Clock drift rate (ms/hr)')")
     *         driftrate
      write(3,"(i5,17x,'Delay in firing guns (integer msec)')")
     *         gundelay
      write(3,"(i5,17x,'Offset(integer msec):obs&ship clocks')")
     *         clockoff
      write(3,"(' SHOT DAY HR MM     SEC   LATITUDE   LONGITUDE',
     *'    RANGE')")
      print *, 'inputing data'
      irec=0
c
   30 read  (2,200,end=99) isht,ida,ihr,min,sec,slong,slat
c  200 format (i4,i4,i3,1x,i2,1x,f6.3,2f11.6)
  200 format (i5,i4,i3,1x,i2,1x,f6.3,2f11.6)
c
	shtlat=slat
	shtlong=slong
      irec=irec+1
      aitime=ida*24+ihr+min/60.+sec/3600.
       if(aitime.lt.stime) goto 30
       if(aitime.gt.etime) goto 99
c
      CALL LLUTM(SHTLAT,SHTLONG,CM,NORTH,EAST)
C      print *, isht,shtlat,shtlong,north,east,shthead
C      pause
      if(shtcor.ne.0.) then
        NORTH=NORTH-SHTCOR*COS(SHTHEAD/RADDEG)
        EAST=EAST-SHTCOR*SIN(SHTHEAD/RADDEG)
      endif
      RNGEAST=EAST-OBSEAST
      RNGNORTH=NORTH-OBSNORTH
      SHTRNG=SQRT(RNGEAST*RNGEAST+RNGNORTH*RNGNORTH)
      SHTRNG=SHTRNG/1000.
	if(rngnorth.lt.0.) shtrng=-shtrng
	if(irec.gt.1) then
	  deln=rngnorth-rngn0
	  dele=rngeast-rnge0
	  if((shtrng-shtrng0).lt.0) shtrng=-shtrng
	  if(deln.gt.dele) then
	    if((rngnorth-rngn0).lt.0) shtrng=-shtrng
	  else
	    if((rngeast-rnge0).lt.0) shtrng=-shtrng
	  endif
      endif
c      CALL UTMLL(SHTLAT,SHTLONG,NORTH,EAST,CM)

      write(3,201) isht,ida,ihr,min,sec,slat,slong,shtrng
  201 format(i5,i4,2i3,f8.4,f11.6,f12.6,f9.3)
      shtrng0=abs(shtrng)
	rnge0=rngeast
	rngn0=rngnorth

      goto 30

   99 continue
      stop
      end
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
