program forward
      implicit double precision (a-h,o-z)
      common/const/pi,rad
      common/elipsoid/a,f
      character*1 answer,option
      character*30 filout,name1,name2,elips
      pi=4.d0*datan(1.d0)
      rad=180.d0/pi
 
    5 write(*,*) '  Program Forward - Version 2.0'
      write(*,*) '  '
      write(*,*) '  Ellipsoid options : '
      write(*,*) '  1) Clark 1866 (NAD27) '
      write(*,*) '  2) GRS80 / WGS84 (NAD83) '
      write(*,*) '  3) Any other ellipsoid '
      write(*,*) '  '
      write(*,*) '  Enter choice : '
      read(*,10) option
   10 format(a1)
      if(option.eq.'1') then
        a=6378206.4d0
        f=1.d0/294.9786982138d0
        elips='Clark 1866     (NAD27)'
      elseif(option.eq.'2') then
        a=6378137.d0
        f=1.d0/298.25722210088d0
        elips='GRS80 / WGS84  (NAD83)'
      elseif(option.eq.'3') then
        call ellips(elips)
      else
        write(*,*) '  Enter 1, 2, or 3 !   Try again --'
        go to 5
      endif
 
   15   write(*,*) '  Enter first station latitude (dd mm ss.sssss)'
        write(*,*) '  Separate d,m,s by blanks or commas : '
        read(*,*) ald1,alm1,sl1
        write(*,*) '  Enter first station longitude (ddd mm ss.sssss)'
        write(*,*) '  Separate d,m,s by blanks or commas : '
        read(*,*) alod1,alom1,slo1
        if(ald1.lt.0.d0) then
          ald1=dabs(ald1)
          alm1=dabs(alm1)
          sl1=dabs(sl1)
          lat1sn=-1
        else
          lat1sn=1
        endif
        if(alod1.lt.0.d0) then
          alod1=dabs(alod1)
          alom1=dabs(alom1)
          slo1=dabs(slo1)
          lon1sn=1
        else
          lon1sn=-1
        endif
        call getrad(ald1,alm1,sl1,lat1sn,glat1)
        call getrad(alod1,alom1,slo1,lon1sn,glon1)
 
   20 write(*,*) '  Enter forward azimuth from north (ddd mm ss.sss)'
      write(*,*) '  Separate d,m,s by blanks or commas : '
      read(*,*) iaz1,maz1,saz1
      iazsn=1
      azd1=dble(iaz1)
      azm1=dble(maz1)
      call getrad(azd1,azm1,saz1,iazsn,faz)
      write(*,*) '  Enter ellipsoidal distance : '
      read(*,*) edist
      call dirct1(glat1,glon1,glat2,glon2,faz,baz,edist)
      call todmsp(faz,iaz1,maz1,saz1,isign1)
      if(isign1.lt.0) then
        iaz1=359-iaz1
        maz1=59-maz1
        saz1=60.d0-saz1
      endif
      call todmsp(baz,iaz2,maz2,saz2,isign2)
      if(isign2.lt.0) then
        iaz2=359-iaz2
        maz2=59-maz2
        saz2=60.d0-saz2
      endif
 
      call todmsp(glat1,ld1,lm1,sl1,lat1sn)
      call todmsp(glon1,lod1,lom1,slo1,lon1sn)
      call todmsp(glat2,ld2,lm2,sl2,lat2sn)
      call todmsp(glon2,lod2,lom2,slo2,lon2sn)
 
      write(*,*) '  '
      write(*,*) '  '
      write(*,*) '  '
      write(*,*) '  '
      write(*,*) '  '
      write(*,49) elips
   49 format('  Ellipsoid : ',a30)
      finv=1.d0/f
      b=a*(1.d0-f)
      write(*,50) a,b,finv
   50 format('  Equatorial radius,  a   = ',f15.4,/,
     *       '  Polar radius,       b   = ',f15.4,/,
     *       '  Inverse flattening, 1/f = ',f16.11)
      write(*,*) '  '
      write(*,*) '  First Station : '
      write(*,*) '  --------------- '
      if(lat1sn.gt.0.d0) then
        write(*,21) ld1,lm1,sl1
   21   format('    LAT = ',i3,1x,i2,1x,f8.5,' North')
      else
        write(*,22) ld1,lm1,sl1
   22   format('    LAT = ',i3,1x,i2,1x,f8.5,' South')
      endif
      if(lon1sn.gt.0.d0) then
        write(*,23) lod1,lom1,slo1
   23   format('    LON = ',i3,1x,i2,1x,f8.5,' East')
      else
        write(*,24) lod1,lom1,slo1
   24   format('    LON = ',i3,1x,i2,1x,f8.5,' West')
      endif
      write(*,*) '  '
      write(*,*) '  Second Station : '
      write(*,*) '  ---------------- '
      if(lat2sn.gt.0.d0) then
        write(*,26) ld2,lm2,sl2
   26   format('    LAT = ',i3,1x,i2,1x,f8.5,' North')
      else
        write(*,27) ld2,lm2,sl2
   27   format('    LAT = ',i3,1x,i2,1x,f8.5,' South')
      endif
      if(lon2sn.gt.0.d0) then
        write(*,28) lod2,lom2,slo2
   28   format('    LON = ',i3,1x,i2,1x,f8.5,' East')
      else
        write(*,29) lod2,lom2,slo2
   29   format('    LON = ',i3,1x,i2,1x,f8.5,' West')
      endif
      write(*,*) '  '
      write(*,32) edist
   32 format('  Ellipsoidal distance  = ',f13.4)
      write(*,34) iaz1,maz1,saz1
   34 format('  Forward azimuth from north = ',i3,1x,i2,1x,f7.4)
      write(*,35) iaz2,maz2,saz2
   35 format('  Back azimuth from north    = ',i3,1x,i2,1x,f7.4)
      write(*,*) '  '
      write(*,*) '  '
      write(*,*) '  '
      write(*,*) '  Do you want to save this output into a file (y/n)?'
      read(*,10) answer
      if(answer.eq.'Y'.or.answer.eq.'y') then
   39   write(*,*) '  Enter the output filename : '
        read(*,40) filout
   40   format(a30)
        open(3,file=filout,status='new',err=99)
        write(3,*) '  '
        write(3,49) elips
        write(3,50) a,b,finv
        write(*,*) '  Enter the first station name : '
        read(*,40) name1
        write(*,*) '  Enter the second station name : '
        read(*,40) name2
        write(3,*) '  '
        write(3,41) name1
   41   format('  First station : ',a30)
        if(lat1sn.gt.0.d0) then
          write(3,21) ld1,lm1,sl1
        else
          write(3,22) ld1,lm1,sl1
        endif
        if(lon1sn.gt.0.d0) then
          write(3,23) lod1,lom1,slo1
        else
          write(3,24) lod1,lom1,slo1
        endif
        write(3,*) '  '
        write(3,42) name2
   42   format('  Second station : ',a30)
        if(lat2sn.gt.0.d0) then
          write(3,26) ld2,lm2,sl2
        else
          write(3,27) ld2,lm2,sl2
        endif
        if(lon2sn.gt.0.d0) then
          write(3,28) lod2,lom2,slo2
        else
          write(3,29) lod2,lom2,slo2
        endif
        write(3,*) '  '
        write(3,32) edist
        write(3,34) iaz1,maz1,saz1
        write(3,35) iaz2,maz2,saz2
        write(3,*) '  '
      endif
      write(*,*) '  1) Another forward, different ellipsoid.'
      write(*,*) '  2) Same ellipsoid, two new stations.'
      write(*,*) '  3) Same ellipsoid, same first station.'
      write(*,*) '  4) Let the second station be the first station.'
      write(*,*) '  5) Done.'
      write(*,*) '  Enter choice : '
      read(*,10) answer
      if(answer.eq.'1') go to 5
      if(answer.eq.'2') go to 15
      if(answer.eq.'3') go to 20
      if(answer.eq.'4') then
        glat1=glat2
        glon1=glon2
        go to 20
      endif
      write(*,*) '  Th Th Th Thats all, folks!'
      stop
   99 write(*,*) '  File already exists, try again.'
      go to 39
      end
 
      subroutine getrad(d,m,s,isign,val)
*** comvert deg, min, sec to radians
      implicit double precision(a-h,j-z)
      common/const/pi,rad
      val=(d+m/60.d0+s/3600.d0)/rad
      val=val*dble(isign)
      return
      end
 
      subroutine todmsp(val,id,im,s,isign)
*** convert position radians to deg,min,sec
*** range is [-pi to +pi]
      implicit double precision(a-h,o-z)
      common/const/pi,rad
    1 if(val.gt.pi) then
        val=val-pi-pi
        go to 1
      endif
    2 if(val.lt.-pi) then
        val=val+pi+pi
        go to 2
      endif
      if(val.lt.0.d0) then
        isign=-1
      else
        isign=+1
      endif
      s=dabs(val*rad)
      id=idint(s)
      s=(s-id)*60.d0
      im=idint(s)
      s=(s-im)*60.d0
*** account for rounding error
      is=idnint(s*1.d5)
      if(is.ge.6000000) then
        s=0.d0
        im=im+1
      endif
      if(im.ge.60) then
        im=0
        id=id+1
      endif
      return
      end
 
      SUBROUTINE DIRCT1(GLAT1,GLON1,GLAT2,GLON2,FAZ,BAZ,S)
C
C *** SOLUTION OF THE GEODETIC DIRECT PROBLEM AFTER T.VINCENTY
C *** MODIFIED RAINSFORD'S METHOD WITH HELMERT'S ELLIPTICAL TERMS
C *** EFFECTIVE IN ANY AZIMUTH AND AT ANY DISTANCE SHORT OF ANTIPODAL
C
C *** A IS THE SEMI-MAJOR AXIS OF THE REFERENCE ELLIPSOID
C *** F IS THE FLATTENING OF THE REFERENCE ELLIPSOID
C *** LATITUDES AND LONGITUDES IN RADIANS POSITIVE NORTH AND EAST
C *** AZIMUTHS IN RADIANS CLOCKWISE FROM NORTH
C *** GEODESIC DISTANCE S ASSUMED IN UNITS OF SEMI-MAJOR AXIS A
C
C *** PROGRAMMED FOR CDC-6600 BY LCDR L.PFEIFER NGS ROCKVILLE MD 20FEB75
C *** MODIFIED FOR SYSTEM 360 BY JOHN G GERGEN NGS ROCKVILLE MD 750608
C
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON/CONST/PI,RAD
      COMMON/ELIPSOID/A,F
      DATA EPS/0.5D-13/
      R=1.-F
      TU=R*DSIN(GLAT1)/DCOS(GLAT1)
      SF=DSIN(FAZ)
      CF=DCOS(FAZ)
      BAZ=0.
      IF(CF.NE.0.) BAZ=DATAN2(TU,CF)*2.
      CU=1./DSQRT(TU*TU+1.)
      SU=TU*CU
      SA=CU*SF
      C2A=-SA*SA+1.
      X=DSQRT((1./R/R-1.)*C2A+1.)+1.
      X=(X-2.)/X
      C=1.-X
      C=(X*X/4.+1)/C
      D=(0.375D0*X*X-1.)*X
      TU=S/R/A/C
      Y=TU
  100 SY=DSIN(Y)
      CY=DCOS(Y)
      CZ=DCOS(BAZ+Y)
      E=CZ*CZ*2.-1.
      C=Y
      X=E*CY
      Y=E+E-1.
      Y=(((SY*SY*4.-3.)*Y*CZ*D/6.+X)*D/4.-CZ)*SY*D+TU
      IF(DABS(Y-C).GT.EPS)GO TO 100
      BAZ=CU*CY*CF-SU*SY
      C=R*DSQRT(SA*SA+BAZ*BAZ)
      D=SU*CY+CU*SY*CF
      GLAT2=DATAN2(D,C)
      C=CU*CY-SU*SY*CF
      X=DATAN2(SY*SF,C)
      C=((-3.*C2A+4.)*F+4.)*C2A*F/16.
      D=((E*CY*C+CZ)*SY*C+Y)*SA
      GLON2=GLON1+X-(1.-C)*D*F
      BAZ=DATAN2(SA,BAZ)+PI
      RETURN
      END

      subroutine ellips(elips)
      implicit double precision(a-h,o-z)
      character*1 answer
      character*30 elips
      common/elipsoid/a,f
      write(*,*) '  Other Ellipsoids.'
      write(*,*) '  -----------------'
      write(*,*) '  '
      write(*,*) '  A) Airy 1858'
      write(*,*) '  B) Airy Modified'
      write(*,*) '  C) Australian National'
      write(*,*) '  D) Bessel 1841'
      write(*,*) '  E) Clark 1880'
      write(*,*) '  F) Everest 1830'
      write(*,*) '  G) Everest Modified'
      write(*,*) '  H) Fisher 1960'
      write(*,*) '  I) Fisher 1968'
      write(*,*) '  J) Hough 1956'
      write(*,*) '  K) International (Hayford)'
      write(*,*) '  L) Krassovsky 1938'
      write(*,*) '  M) NWL-9D (WGS 66)'
      write(*,*) '  N) South American 1969'
      write(*,*) '  O) Soviet Geod. System 1985'
      write(*,*) '  P) WGS 72'
      write(*,*) '  Q-Z) You define it.'
      write(*,*) '  '
      write(*,*) '  Enter choice : '
      read(*,10) answer
   10 format(a1)
      if(answer.eq.'A'.or.answer.eq.'a') then
        a=6377563.396d0
        f=1.d0/299.3249646d0
        elips='Airy 1858'
      elseif(answer.eq.'B'.or.answer.eq.'b') then
        a=6377340.189d0
        f=1.d0/299.3249646d0
        elips='Airy Modified'
      elseif(answer.eq.'C'.or.answer.eq.'c') then
        a=6378160.d0
        f=1.d0/298.25d0
        elips='Australian National'
      elseif(answer.eq.'D'.or.answer.eq.'d') then
        a=6377397.155d0
        f=1.d0/299.1528128d0
        elips='Bessel 1841'
      elseif(answer.eq.'E'.or.answer.eq.'e') then
        a=6378249.145d0
        f=1.d0/293.465d0
        elips='Clarke 1880'
      elseif(answer.eq.'F'.or.answer.eq.'f') then
        a=6377276.345d0
        f=1.d0/300.8017d0
        elips='Everest 1830'
      elseif(answer.eq.'G'.or.answer.eq.'g') then
        a=6377304.063d0
        f=1.d0/300.8017d0
        elips='Everest Modified'
      elseif(answer.eq.'H'.or.answer.eq.'h') then
        a=6378166.d0
        f=1.d0/298.3d0
        elips='Fisher 1960'
      elseif(answer.eq.'I'.or.answer.eq.'i') then
        a=6378150.d0
        f=1.d0/298.3d0
        elips='Fisher 1968'
      elseif(answer.eq.'J'.or.answer.eq.'j') then
        a=6378270.d0
        f=1.d0/297.d0
        elips='Hough 1956'
      elseif(answer.eq.'K'.or.answer.eq.'k') then
        a=6378388.d0
        f=1.d0/297.d0
        elips='International (Hayford)'
      elseif(answer.eq.'L'.or.answer.eq.'l') then
        a=6378245.d0
        f=1.d0/298.3d0
        elips='Krassovsky 1938'
      elseif(answer.eq.'M'.or.answer.eq.'m') then
        a=6378145.d0
        f=1.d0/298.25d0
        elips='NWL-9D  (WGS 66)'
      elseif(answer.eq.'N'.or.answer.eq.'n') then
        a=6378160.d0
        f=1.d0/298.25d0
        elips='South American 1969'
      elseif(answer.eq.'O'.or.answer.eq.'o') then
        a=6378136.d0
        f=1.d0/298.257d0
        elips='Soviet Geod. System 1985'
      elseif(answer.eq.'P'.or.answer.eq.'p') then
        a=6378135.d0
        f=1.d0/298.26d0
        elips='WGS 72'
      else
        write(*,*) '  Enter equatorial radius, a : '
        read(*,*) a
        b=0.d0
        write(*,*) '  Enter polar radius, b.  If you want to enter'
        write(*,*) '  flattening instead, enter 0 (zero) : '
        read(*,*) b
        if(b.lt.0.1d0) then
          write(*,*) '  Enter 1/f : '
          read(*,*) finv
          f=1.d0/finv
        else
          f=(a-b)/a
        endif
        elips='User defined.'
      endif
      return
      end