PROGRAM FITPAR
C
      REAL*8 AZMTH(100),HGHT(100),Z0
      REAL*8 SUMX,SUMX2,SUMX3,SUMX4,SUMXY,SUMY,SUMX2Y
      REAL*8 DET1,DET2,DET3,DET,A,B,C,XXX
      REAL*8 TEMP,PRESSURE,REFR
      CHARACTER*80 INFILE,OUTFILE
C
      DTAND(XXX)=DTAN(XXX*3.1415926D0/180.0D0)
C
C     Pressure in milibars
      PRESSURE=1013.25D0
c     Temperature in Celsius
      TEMP=20.0D0
C
c     input data
      PRINT *, 'ENTER INPUT FILE W/ H AND A VALUES:'
      READ(*,'(A)') INFILE
      PRINT *, 'ENTER OUTPUT FILE :'
      READ(*,'(A)') OUTFILE
      PRINT *, 'ENTER ZENITAL ERROR (IN ARCMIN)'
      READ(*,*) Z0
      Z0=Z0/60.D0
c
c     choose input format options 
      PRINT *, 'TYPE 1 FOR D:M:S INSTEAD OF DECIMAL DEG:'
      READ(*,*) IOP0
C
      PRINT *, 'TYPE 1 FOR ZENITAL DIST. INSTEAD OF ALTITUDE:'
      READ(*,*) IOP1
C
      PRINT *, 'TYPE 1 FOR (H,V) ORDER INSTEAD OF (V,H):'
      READ(*,*) IOP2
C
      OPEN (1,FILE=INFILE)
      OPEN (2,FILE=OUTFILE)
C
      SUMX=0.
      SUMX2=0.
      SUMX3=0.
      SUMX4=0.
      SUMXY=0.
      SUMY=0.
      SUMX2Y=0.
      DO I=1,100
         IF (IOP0.NE.1) THEN
            IF (IOP2.NE.1) THEN
               READ(1,*,END=10) HGHT(I),AZMTH(I)
            ELSE
               READ(1,*,END=10) AZMTH(I),HGHT(I)
            ENDIF
         ELSE
            IF (IOP2.NE.1) THEN
               READ(1,*,END=10) HD,HM,HS,AD,AM,AS
            ELSE
               READ(1,*,END=10) AD,AM,AS,HD,HM,HS
            ENDIF
            IF (HD.GE.0.) THEN
               HGHT(I) = HD + HM/60. + HS/3600.
            ELSE
               HGHT(I) = HD - HM/60. - HS/3600.
            ENDIF
            IF (AD.GE.0.) THEN
               AZMTH(I) = AD + AM/60. + AS/3600.
            ELSE
               AZMTH(I) = AD - AM/60. - AS/3600.
            ENDIF
         ENDIF            
C     CORRECT Z FOR ZENITAL ERROR
         IF (IOP1.EQ.1) THEN
            HGHT(I) = HGHT(I) - Z0 
            HGHT(I) = 90.D0 - HGHT(I)
         ELSE
            HGHT(I) = HGHT(I) + Z0
         ENDIF 
C     CORRECT FOR ATMOSPHERIC REFRACTION
         REFR = 0.00452D0*PRESSURE/((TEMP + 273.0D0)*DTAND(HGHT(I)))
         HGHT(I) = HGHT(I) - REFR
C
         WRITE(2,*) AZMTH(I),HGHT(I)
C
C     AZMTH(I)=AZMTH(I)/100.
C     HGHT(I)=HGHT(I)/100.
         SUMX=SUMX+AZMTH(I)
         SUMX2=SUMX2+AZMTH(I)*AZMTH(I)
         SUMX3=SUMX3+AZMTH(I)*AZMTH(I)*AZMTH(I)
         SUMX4=SUMX4+AZMTH(I)*AZMTH(I)*AZMTH(I)*AZMTH(I)
         SUMY=SUMY+HGHT(I)
         SUMXY=SUMXY+AZMTH(I)*HGHT(I)
         SUMX2Y=SUMX2Y+AZMTH(I)*AZMTH(I)*HGHT(I)
c         print *, 'az,h:',azmth(i),hght(i)
c         print *, 'sumx,sumx2,sumx3,sumx4:',sumx,sumx2,sumx3,sumx4
c         print *, 'sumy,sumxy,sumx2y:',sumy,sumxy,sumx2y
c         pause 'check:'
      END DO
 10   CLOSE (1)
      CLOSE (2)
      NPTS=I-1
      PRINT *, 'NUMBER OF DATA POINTS:',NPTS
      DET=NPTS*SUMX2*SUMX4 + SUMX*SUMX3*SUMX2 + SUMX2*SUMX*SUMX3 - 
     &     SUMX2*SUMX2*SUMX2 - SUMX3*SUMX3*NPTS - SUMX4*SUMX*SUMX
      DET1=SUMY*SUMX2*SUMX4 + SUMX*SUMX3*SUMX2Y + SUMX2*SUMXY*SUMX3 - 
     &     SUMX2Y*SUMX2*SUMX2 - SUMX3*SUMX3*SUMY - SUMX4*SUMXY*SUMX
      DET2 = NPTS*SUMXY*SUMX4 + SUMY*SUMX3*SUMX2 + SUMX2*SUMX*SUMX2Y -
     &     SUMX2*SUMXY*SUMX2 - SUMX2Y*SUMX3*NPTS - SUMX4*SUMX*SUMY
      DET3=NPTS*SUMX2*SUMX2Y + SUMX*SUMXY*SUMX2 + SUMY*SUMX*SUMX3 - 
     &     SUMX2*SUMX2*SUMY - SUMX3*SUMXY*NPTS - SUMX2Y*SUMX*SUMX
      PRINT *, 'DET,DET1,DET2,DET3:',DET,DET1,DET2,DET3
      C = DET1/DET
      B = DET2/DET
      A = DET3/DET
C
      PRINT *, 'BEST FIT PARS:',A,B,C
C      
      STOP
      END