      PROGRAM VRADIAL
C CALCULATES AND DISPLAYS RADIAL VELOCITY DATA AND MODELS
C
C  3 Sept 1991        J.T. Armstrong

      IMPLICIT UNDEFINED (A-Z)

      INCLUDE 'VELPARM.INC'

      REAL*8 A1_SINI,A2_SINI,V0,DV0DT
      REAL*8 TOBS(MAXVDATA),V1OBS(MAXVDATA),V2OBS(MAXVDATA)
      REAL*8 V1MOD(MAXVDATA),V2MOD(MAXVDATA)
      REAL*8 DELV1(MAXVDATA),DELV2(MAXVDATA),DELVSQ
      REAL*8 CHISQ,CHISQ1,CHISQ2
      REAL*8 TCURV(MAXVMOD),V1CURV(MAXVMOD),V2CURV(MAXVMOD)
      REAL*8 T,A(11),V1,V2,DVELDA(11,2)
      REAL*8 E_TRUE,PERIOD,EPOCH,ARG_PERI
      REAL*8 DV1DA1,DV2DA2,DV1DV0,DV2DV0
      REAL*8 NEWVALUE
      REAL*8 SIG1(MAXVDATA),SIG2(MAXVDATA)
      REAL*8 PUNCERT3D(MMAX,-1:1)

      INTEGER*4 I,J,FK5,INDEX,NVPTS,IERR,BADPTS,BADPTS1,BADPTS2
      INTEGER*4 GOODPTS, GOODPTS1, GOODPTS2

      CHARACTER*1 ANSWER
      CHARACTER*128 VELFILE, OUTFILE
      CHARACTER*80 STARNAME, COMMENT
      CHARACTER*28 ASTRING

      LOGICAL FIRST, DO_1, DO_2

      FIRST = .TRUE.
      DO_1 = .TRUE.
      DO_2 = .TRUE.

C GET DATA FILE NAME AND READ IT

      WRITE(OUTC,'(A)') 
     1    ' VRADIAL takes file w/ V1, sig(V1), V2, sig(V2), and fits'
      WRITE(OUTC,'(A)') 
     1    ' a spectroscopic orbit using a1 sin i, a2 sin i, e, P, T,'
      WRITE(OUTC,'(A)') 
     1    ' omega, V(0), and dV(0)/dt (in case there is a third body).'
      WRITE(OUTC,'(A)') ' ENTER DATA FILE NAME: '
      READ(INC,*) VELFILE
      OPEN(UNIT=INVDATA, FILE=VELFILE, FORM='FORMATTED',STATUS='OLD')
      CALL GETVEL(INDEX,STARNAME,FK5,PERIOD,TOBS,V1OBS,SIG1,
     1            V2OBS,SIG2,INVDATA)
      CLOSE (UNIT=INVDATA)

      NVPTS = INDEX - 2
      WRITE(OUTC,'(A,I4)') ' VRADIAL:  NVPTS = ',NVPTS
      WRITE(OUTC,'(1X,A,A30,A,I4,A,F8.3)')
     1      'STAR ',STARNAME,'; FK5',FK5,'; PERIOD ',PERIOD
      IF(NVPTS.GT.MAXVDATA) THEN
          WRITE(OUTC,'(A,I3,A,I3)') ' NVPTS = ',NVPTS,
     1                ' EXCEEDS MAXVDATA = ',MAXVDATA
          GO TO 999                   ! GO TO EXIT
      END IF

C COUNT UP THE BAD POINTS (V = 999)
      BADPTS = 0
      BADPTS1 = 0
      BADPTS2 = 0
      GOODPTS = 0
      GOODPTS1 = 0
      GOODPTS2 = 0
      DO I=1,NVPTS
          IF(V1OBS(I).GE.999.) THEN
              BADPTS=BADPTS+1
              BADPTS1=BADPTS1+1
          ELSE
              GOODPTS=GOODPTS+1
              GOODPTS1=GOODPTS1+1
          END IF
          IF(V2OBS(I).GE.999.) THEN
              BADPTS=BADPTS+1
              BADPTS2=BADPTS2+1
          ELSE
              GOODPTS=GOODPTS+1
              GOODPTS2=GOODPTS2+1
          END IF
      END DO

C ASK FOR INPUT PARAMETERS
C RETURN HERE IF ANOTHER MODEL PLOT IS WANTED

  100 IF(FIRST) THEN
          WRITE(OUTC,'(A)') 
     1    ' ENTER a1*sini, a2*sini (KM), e, P, T, omega, V0, dV0/dt: '
          READ(INC,*) A1_SINI, A2_SINI, E_TRUE, PERIOD, EPOCH,
     1                ARG_PERI, V0, DV0DT
          ARG_PERI = ARG_PERI*RPDEG
  101     WRITE(OUTC,*) ' Do V1, V2, or both? [1,2,B]:'
          READ(INC,*) ANSWER
          IF(ANSWER.EQ.'1') THEN
              DO_1 = .TRUE.
              DO_2 = .FALSE.
          ELSE IF(ANSWER.EQ.'2') THEN
              DO_1 = .FALSE.
              DO_2 = .TRUE.
          ELSE IF(ANSWER.EQ.'B'.OR.ANSWER.EQ.'b') THEN
              DO_1 = .TRUE.
              DO_2 = .TRUE.
          ELSE
              GO TO 101
          END IF
          FIRST = .FALSE.
          A(1) = 0.
          A(2) = E_TRUE
          A(3) = PERIOD
          A(4) = EPOCH
          A(5) = 0.
          A(6) = 0.
          A(7) = ARG_PERI
          A(8) = A1_SINI
          A(9) = A2_SINI
          A(10)= V0
          A(11)= DV0DT
      ELSE
          WRITE(OUTC,*) ' Enter command: '
          WRITE(OUTC,'(A,A)') 
     1      ' Change parameter: 1, 2 [a1,2 sin i] e, P, T, w,',
     2      ' V, A(cceleration), N [none]'
          WRITE(OUTC,'(A)') ' C [change stars])'
          WRITE(OUTC,'(A)') ' F [fit velocity curve]'
          WRITE(OUTC,'(A)') ' R [radial velocity for an input date]'
          READ(INC,'(A)') ANSWER
          CALL CAPS(ANSWER)
          IF(ANSWER.EQ.'1' .OR. ANSWER.EQ.'2' .OR. ANSWER.EQ.'E'
     1       .OR. ANSWER.EQ.'P' .OR. ANSWER.EQ.'T' .OR. ANSWER.EQ.'W'
     2       .OR. ANSWER.EQ.'V' .OR. ANSWER.EQ.'A') THEN
              WRITE(OUTC,*) ' Enter new value for ',ANSWER,': '
              READ(INC,*) NEWVALUE
              IF(ANSWER.EQ.'1') A1_SINI = NEWVALUE
              IF(ANSWER.EQ.'2') A2_SINI = NEWVALUE
              IF(ANSWER.EQ.'E') E_TRUE  = NEWVALUE
              IF(ANSWER.EQ.'P') PERIOD  = NEWVALUE
              IF(ANSWER.EQ.'T') EPOCH   = NEWVALUE
              IF(ANSWER.EQ.'W') ARG_PERI= NEWVALUE*RPDEG
              IF(ANSWER.EQ.'V') V0      = NEWVALUE
              IF(ANSWER.EQ.'A') DV0DT   = NEWVALUE
              GO TO 100                       ! GET ANOTHER NEW VALUE
          ELSE IF(ANSWER.EQ.'C') THEN
  102         WRITE(OUTC,*) ' Do V1, V2, or both? [1,2,B]:'
              READ(INC,*) ANSWER
              IF(ANSWER.EQ.'1') THEN
                  DO_1 = .TRUE.
                  DO_2 = .FALSE.
              ELSE IF(ANSWER.EQ.'2') THEN
                  DO_1 = .FALSE.
                  DO_2 = .TRUE.
              ELSE IF(ANSWER.EQ.'B'.OR.ANSWER.EQ.'b') THEN
                  DO_1 = .TRUE.
                  DO_2 = .TRUE.
              ELSE
                  GO TO 102
              END IF
              GO TO 100                       ! GET NEW VALUES
          ELSE IF(ANSWER.EQ.'R') THEN
              WRITE(OUTC,*) ' Type JD - 2440000:'
              READ(INC,*) NEWVALUE
              NEWVALUE = NEWVALUE + 2440000.0
          END IF                       

          A(1) = 0.
          A(2) = E_TRUE
          A(3) = PERIOD
          A(4) = EPOCH
          A(5) = 0.
          A(6) = 0.
          A(7) = ARG_PERI
          A(8) = A1_SINI
          A(9) = A2_SINI
          A(10)= V0
          A(11)= DV0DT

          IF(ANSWER.EQ.'F') THEN
              CALL MRQSHELL(TOBS,V1OBS,V2OBS,SIG1,SIG2,NVPTS,A,DELVSQ,
     1                      STARNAME,DO_1,DO_2,PUNCERT3D,BADPTS)
              E_TRUE  = A(2) 
              PERIOD  = A(3) 
              EPOCH   = A(4) 
              ARG_PERI= A(7) 
              A1_SINI = A(8) 
              A2_SINI = A(9) 
              V0      = A(10)
              DV0DT   = A(11)
          ELSE IF(ANSWER.EQ.'R') THEN
              CALL TRUE2VEL(NEWVALUE,A,V1,V2,DVELDA)
              WRITE(OUTC,*)' For JD           V1        V2'
              WRITE(OUTC,'(F12.3,F10.3,F10.3)') NEWVALUE,V1,V2
              GO TO 100
          END IF

      END IF

C CALCULATE MODEL VELOCITY POINTS FOR COMPARISON TO DATA
C RETURN HERE TO REPEAT THE MODEL

  200 DELVSQ=0.
      CHISQ=0.
      CHISQ1=0.
      CHISQ2=0.

      DO I=1,NVPTS
          T=TOBS(I)
          CALL TRUE2VEL(T,A,V1,V2,DVELDA)
          V1MOD(I)=V1
          V2MOD(I)=V2
          DELV1(I)=V1OBS(I)-V1MOD(I)
          DELV2(I)=V2OBS(I)-V2MOD(I)
          IF(DO_1 .AND. (V1OBS(I).LT.999.) ) THEN
              DELVSQ = DELVSQ + DELV1(I)**2
              CHISQ=CHISQ+(DELV1(I)/SIG1(I))**2
              CHISQ1=CHISQ1+(DELV1(I)/SIG1(I))**2
          END IF
          IF(DO_2 .AND. (V2OBS(I).LT.999.) ) THEN
              DELVSQ = DELVSQ + DELV2(I)**2
              CHISQ=CHISQ+(DELV2(I)/SIG2(I))**2
              CHISQ2=CHISQ2+(DELV2(I)/SIG2(I))**2
          END IF
      END DO
      CHISQ=CHISQ/GOODPTS
      CHISQ1=CHISQ1/GOODPTS1
      CHISQ2=CHISQ2/GOODPTS2

C CALCULATE MODEL VELOCITY POINTS FOR PLOTTING ON THE SCREEN

      DO I=1,MAXVMOD
          TCURV(I)=EPOCH+PERIOD*(I-1)/(MAXVMOD-1)
          T=TCURV(I)
          CALL TRUE2VEL(T,A,V1,V2,DVELDA)
          V1CURV(I)=V1
          V2CURV(I)=V2
      END DO

C PLOT RESULT ON SCREEN

      CALL PLOTVEL(A,TOBS,TCURV,V1OBS,V1CURV,V2OBS,V2CURV,V1MOD,
     1             V2MOD,DELV1,DELV2,CHISQ,
     2             NVPTS,STARNAME,FK5,DO_1,DO_2)

C LIST RESULT TO SCREEN

      IF(DO_1.AND..NOT.DO_2) ASTRING='ONLY STAR 1 DONE;   CHI^2 = '
      IF(DO_2.AND..NOT.DO_1) ASTRING='ONLY STAR 2 DONE;   CHI^2 = '
      IF(DO_1.AND.DO_2)      ASTRING='STAR 1 AND STAR 2:  CHI^2 = '

      WRITE(OUTC,*) ' List result to screen?'
      READ(INC,*) ANSWER
      CALL CAPS(ANSWER)
      IF(ANSWER.EQ.'Y') THEN
          WRITE(OUTC,'(A,A)')
     1     '  a1*sin i   a2*sin i   eccen.    period',
     2     '       epoch     arg peri   V0   dV0/dt'
          WRITE(OUTC,'(1X,E10.3,1X,E10.3,1X,F8.5,1X,F10.5,1X,F11.3,1X,
     1           F8.3,1X,F7.3,1X,F7.4)')
     2           A1_SINI,A2_SINI,E_TRUE,PERIOD,EPOCH,ARG_PERI/RPDEG,V0,
     3           DV0DT
          WRITE(OUTC,'(1X,A,F10.3)') ASTRING,CHISQ 
          IF(DO_1) WRITE(OUTC,*) ' STAR 1 CHI^2 = ', CHISQ1
          IF(DO_2) WRITE(OUTC,*) ' STAR 2 CHI^2 = ', CHISQ2

C RULER LINE:
CaaaaaaaaaaXaaaaaaaaaaXeeeeeeeeXppppppppppXtttttttttttttXooooooooXvvvvvvv

          WRITE(OUTC,'(A)')
     1  '   JD obs     v1 obs  v2 obs  v1 mod  v2 mod   v1 (O-C) v2'

C RULER LINE:
CXtttttttttttXXv1v1v1XXv2v2v2XXv1v1v1XXv2v2v2XXv1v1v1XXv2v2v2

          DO I=1,NVPTS
              WRITE(OUTC,'(1X,F11.3,6(2X,F6.2))')
     1            TOBS(I),V1OBS(I),V2OBS(I),V1MOD(I),
     2            V2MOD(I),DELV1(I),DELV2(I)
          END DO
      END IF

C WRITE RESULTS TO OUTPUT FILE:

      WRITE(OUTC,*) 'Write result to file?' 
      READ(INC,*) ANSWER
      CALL CAPS(ANSWER)
      IF(ANSWER.EQ.'Y') THEN
  300     WRITE(OUTC,'(A)') ' ENTER OUTPUT FILE NAME:'
          READ(INC,*) OUTFILE
          OPEN(UNIT=OUTVFILE, FILE=OUTFILE, FORM='FORMATTED',
     1         STATUS='NEW',IOSTAT=IERR)
          IF(IERR.NE.0) THEN
              WRITE(OUTC,*) ' Cannot open file ',OUTFILE
              GO TO 300
          END IF
          WRITE(OUTC,*) 'ENTER COMMENT:'
          READ(INC,'(A80)') COMMENT
          WRITE(OUTVFILE,'(1X,A,F10.3)') ASTRING,CHISQ
          IF(DO_1) WRITE(OUTC,*) ' STAR 1 CHI^2 = ', CHISQ1
          IF(DO_2) WRITE(OUTC,*) ' STAR 2 CHI^2 = ', CHISQ2
          WRITE(OUTVFILE,'(A)') COMMENT
          WRITE(OUTVFILE,'(A,A)')
     1         '  a1*sin i   a2*sin i   eccen.    period',
     2         '       epoch     arg peri   V0  dV0/dt'
          WRITE(OUTVFILE,'(1X,E10.3,1X,E10.3,1X,F8.5,1X,F10.5,1X,F11.3,
     1           1X,F8.3,1X,F7.3,1X,F7.4)')
     2           A1_SINI,A2_SINI,E_TRUE,PERIOD,EPOCH,ARG_PERI/RPDEG,V0,
     3           DV0DT
          WRITE(OUTVFILE,'(A)')
     1  '   JD obs     v1 obs  v2 obs  v1 mod  v2 mod   v1 (O-C) v2'
          DO I=1,NVPTS
              WRITE(OUTVFILE,'(1X,F11.3,6(2X,F6.2))')
     1        TOBS(I),V1OBS(I),V2OBS(I),V1MOD(I),
     2        V2MOD(I),DELV1(I),DELV2(I)
          END DO
          CLOSE(OUTVFILE)
      END IF

C TRY AGAIN?

      WRITE(OUTC,*) ' TRY ANOTHER SET OF PARAMETERS?'
      WRITE(OUTC,*) ' [Y = yes; R = repeat; N = no]'
      READ(INC,*) ANSWER
      CALL CAPS(ANSWER)
      IF(ANSWER.EQ.'Y') GO TO 100
      IF(ANSWER.EQ.'R') GO TO 200

C COME HERE IF NVPTS > MAXVDATA

  999 CONTINUE

      STOP
      END


