PROGRAM GPS_2_DIG

    ! A simple utility program to allow plotting of geodetic benchmark
    ! locations (or any other data points beginning with {lon, lat})
    ! in interactive F-E grid editors OrbWeaver and/or OrbWin,
    ! by reading a *.GPS file (Peter Bird's personal format),
    ! and producing a *.DIG file (Peter Bird's personal format)
    ! which shows the benchmark location with a small cross (+) and triangle.
    ! By Peter Bird, UCLA, 2002.09.06; updated 2013.01.03 & 2015.07.12

    USE DSphere ! spherical trigonometry, lon, lat, and uvec functions;
                ! Fortran 90 MODULE DSphere, in Peter Bird's file DSphere.f90

    IMPLICIT NONE
    CHARACTER*80    :: gps_file, dig_file
    CHARACTER*132   :: gps_format, line
    INTEGER         :: ios
    REAL*8, PARAMETER :: R_in_km = 6371.0D0
    REAL*8            :: lon, lat, &
                       & r_radians, &
                       & symbol_in_km, symbol_in_radians, &
                       & t_radians
    REAL*8, DIMENSION(3) :: omega_uvec, uvec, uvec1, uvec2, uvec3

    WRITE (*, "( ' ')")
    WRITE (*, "( ' ')")

    WRITE (*, "(/' PROGRAM GPS_2_DIG')")
    WRITE (*, "(/' A simple utility program to allow plotting of geodetic benchmarks')")
    WRITE (*, "( ' (or any other data points beginning with {lon, lat})')")
    WRITE (*, "( ' in interactive F-E grid editors OrbWeaver and/or OrbWin,')")
    WRITE (*, "( ' by reading a *.GPS file (Peter Bird''s personal format),')")
    WRITE (*, "( ' and producing a *.DIG file (Peter Bird''s personal format)')")
    WRITE (*, "( ' which shows the benchmark location with a small cross (+) and triangle.')")
    
 1  WRITE (*, "(/' Enter name of existing .GPS (or other type of point-data) file: ')")
    READ  (*, "(A)") gps_file
    OPEN (UNIT = 1, FILE = gps_file, STATUS = 'OLD', IOSTAT = ios)
    IF (ios /= 0) THEN
        WRITE (*, "(' ERROR: File not found (in current directory).')")
        CALL Pause()
        GO TO 1
    END IF

    WRITE (*, *)
    WRITE (*, "(' Caution: If data file is not in .GPS format, first 3 lines will be ignored!')")
    
 2  WRITE (*, "(/' Enter name for new .DIG file: ')")
    READ  (*, "(A)") dig_file
    OPEN (UNIT = 2, FILE = dig_file, STATUS = 'NEW', IOSTAT = ios)
    IF (ios /= 0) THEN
        WRITE (*, "(' ERROR: File already exists (in current directory).')")
        CALL Pause()
        GO TO 2
    END IF

    WRITE (*, "(/' Enter desired width of triangle symbol, in kilometers (e.g., 4): '\)")
    READ (*, *) symbol_in_km
    symbol_in_radians = symbol_in_km / R_in_km
    r_radians = 0.57735D0 * symbol_in_radians ! size of radial step in forming triangle
    t_radians = r_radians / 3.0D0 ! size of sidestep in forming +

    READ (1, "(A)") line ! file name and comments on source
    READ (1, "(A)") GPS_format
    READ (1, *) ! throw away column headers
    counting: DO 
        READ (1, *, IOSTAT = ios) lon, lat ! rest of record ignored (may NOT be .gps format!)
        IF (ios /= 0) EXIT counting

        ! center point:
        CALL DLonLat_2_Uvec (lon, lat, uvec)

        ! triangle symbol:
        CALL DTurn_To (azimuth_radians = 0.0D0, &
                     & base_uvec = uvec, &
                     & far_radians = r_radians, & ! inputs
                     & omega_uvec = omega_uvec, result_uvec = uvec1) ! outputs
        CALL DUvec_2_LonLat (uvec1, lon, lat)
 20     FORMAT (1X,SP,ES12.5,',',ES12.5)
        WRITE (2, 20) lon, lat

        CALL DTurn_To (azimuth_radians = 2.0944D0, &
                     & base_uvec = uvec, &
                     & far_radians = r_radians, & ! inputs
                     & omega_uvec = omega_uvec, result_uvec = uvec2) ! outputs
        CALL DUvec_2_LonLat (uvec2, lon, lat)
        WRITE (2, 20) lon, lat

        CALL DTurn_To (azimuth_radians = 4.189D0, &
                    & base_uvec = uvec, &
                    & far_radians = r_radians, & ! inputs
                    & omega_uvec = omega_uvec, result_uvec = uvec3) ! outputs
        CALL DUvec_2_LonLat (uvec3, lon, lat)
        WRITE (2, 20) lon, lat

        CALL DUvec_2_LonLat (uvec1, lon, lat)
        WRITE (2, 20) lon, lat

  21    FORMAT ('*** end of line segment ***')
        WRITE (2, 21)

        ! + symbol:
        CALL DTurn_To (azimuth_radians = 0.0D0, &
                     & base_uvec = uvec, &
                     & far_radians = t_radians, & ! inputs
                     & omega_uvec = omega_uvec, result_uvec = uvec1) ! outputs
        CALL DUvec_2_LonLat (uvec1, lon, lat)
        WRITE (2, 20) lon, lat

        CALL DTurn_To (azimuth_radians = 3.14159D0, &
                     & base_uvec = uvec, &
                     & far_radians = t_radians, & ! inputs
                     & omega_uvec = omega_uvec, result_uvec = uvec2) ! outputs
        CALL DUvec_2_LonLat (uvec2, lon, lat)
        WRITE (2, 20) lon, lat
        WRITE (2, 21)
        CALL DTurn_To (azimuth_radians = 1.5708D0, &
                     & base_uvec = uvec, &
                     & far_radians = t_radians, & ! inputs
                     & omega_uvec = omega_uvec, result_uvec = uvec1) ! outputs
        CALL DUvec_2_LonLat (uvec1, lon, lat)
        WRITE (2, 20) lon, lat

        CALL DTurn_To (azimuth_radians = 4.7124D0, &
                     & base_uvec = uvec, &
                     & far_radians = t_radians, & ! inputs
                     & omega_uvec = omega_uvec, result_uvec = uvec2) ! outputs
        CALL DUvec_2_LonLat (uvec2, lon, lat)
        WRITE (2, 20) lon, lat
        WRITE (2, 21)

    END DO counting

    CLOSE (1)
    CLOSE (2)
    WRITE (*, "(/' Job completed.')")
    WRITE (*, "(/' NOTE that if you wish to see benchmarks in relation to fault traces in Orbwin,')")
    WRITE (*, "( ' you can simply concatenate this new .DIG file with any existing fault-trace')")
    WRITE (*, "( ' f_*.DIG file, to create a temporary merged .DIG file for viewing in OrbWin.')")
    CALL Pause()

CONTAINS

    SUBROUTINE Pause()
       IMPLICIT NONE
       WRITE (*,"(' Press [Enter]...'\)")
       READ  (*,*)
    END SUBROUTINE Pause

END PROGRAM gps_2_dig