PROGRAM GPS_2_DIG

    ! A simple utility program to plot geodetic benchmarks in OrbWeaver,
    ! by reading a *.gpgs 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

    USE Sphere ! spherical trigonometry, lon, lat, and uvec functions

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

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

    WRITE (*, "(/' PROGRAM GPS_2_DIG')")
    WRITE (*, "(/' A simple utility program to plot geodetic benchmarks in OrbWeaver,')")
    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 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
    
 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.57735 * symbol_in_radians ! size of radial step in forming triangle
    t_radians = r_radians / 3.0 ! 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, gps_format, IOSTAT = ios) &
             & lon, lat, &
             & benchmark_E_velocity, benchmark_N_velocity, &
             & benchmark_E_sigma, benchmark_N_sigma, &
             & benchmark_correlation
        IF (ios /= 0) EXIT counting

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

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

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

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

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

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

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

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

        CALL Turn_To (azimuth_radians = 4.7124, &
                    & base_uvec = uvec, &
                    & far_radians = t_radians, & ! inputs
                    & omega_uvec = omega_uvec, result_uvec = uvec2) ! outputs
        CALL Uvec_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