' CALIBRAT.BAS - TURBO BASIC program to calibrate an EDM instrument run over
' a base line, using equations developed in NOAA Technical Memorandum
' NOS NGS-10 pages 10-23.
'
' Both the Published Base Line Distances File and User Observations File
' must have records of the form:
'          FROM point number, TO point number, Distance
' with the numbers seperated by either commas or blanks.
' The records can be in any order.
' The point numbers need only be unique; they do not need to be consecutive.
' (i.e. You can use point designations for the point number.)
'
' Constraints: 2 < Number of Base Line Points < 16
'              2 < Number of User Observations < 211
'
' A null response for either the Published Base Line Distances Filename or
' User Observations Filename terminates execution (way to bail out).
'
' The results are written to file RESULTS.DAT


DEFDBL A-Z
CLS

OPEN "RESULTS.DAT" FOR OUTPUT AS #2

DIM Delta(210)

INPUT "How many points are on this baseline? : " , Response
NumberOfPoints% = FIX(Response)
WHILE Response <> NumberOfPoints% OR NumberOfPoints% < 3 OR NumberOfPoints% >15
  PRINT
  IF Response <> NumberOfPoints% THEN
    PRINT "I have a problem with that number.  Try again!"
  ELSEIF NumberOfPoints% < 3 THEN
    PRINT "The base line MUST have at least 3 points.  Try again!"
  ELSE
    PRINT "The base line CAN NOT have more than 15 points.  ";_
          "Try again!"
  END IF
  PRINT:PRINT
  INPUT "How many points are on this baseline? : " , Response
  NumberOfPoints% = FIX(Response)
WEND
BaseSize% = (NumberOfPoints% * (NumberOfPoints% - 1)) / 2

DIM Published(1:BaseSize%,1:3)

ON ERROR GOTO BadFile
DO
  NoGood$ = "F"
  PRINT:PRINT
  PRINT "What is the fully qualified name of the file"
  INPUT "   containing the Published Distances? : ", PUBDAT$
  IF PUBDAT$ = "" THEN
    PRINT
    PRINT"I will take a null file name as a signal to stop executing!"
    CLOSE #2
    END
  END IF
  OPEN PUBDAT$ FOR INPUT AS #1
LOOP WHILE NoGood$ = "T"

ON ERROR GOTO PubEOF
WHILE NOT EOF(1)
  X1 = 0.0
  X2 = 0.0
  X3 = 0.0
  INPUT #1, X1, X2, X3
  IF X1 <> FIX(X1) THEN
    PRINT
    PRINT "The following record has a bad FROM point number"
    PRINT "in the Published Distances File:"
    PRINT X1, X2, X3
    PRINT #2," "
    PRINT #2,"The following record has a bad FROM point number"
    PRINT #2,"in the Published Distances File:"
    PRINT #2, X1 ","  X2 "," X3
  ELSEIF X2 <> FIX(X2) THEN
    PRINT
    PRINT "The following record has a bad TO point number"
    PRINT "in the Published Distances File:"
    PRINT X1, X2, X3
    PRINT #2," "
    PRINT #2,"The following record has a bad TO point number"
    PRINT #2,"in the Published Distances File:"
    PRINT #2, X1 "," X2 "," X3
  ELSE
    I% =I% + 1
    IF I% <= BaseSize% THEN
      Published(I%,1) = X1
      Published(I%,2) = X2
      Published(I%,3) = X3
    END IF
  END IF
WEND
CLOSE #1
ON ERROR GOTO 0

IF I% <> BaseSize% THEN
  PRINT
  PRINT USING "For a ## point baseline there must be ### published distances";_
               NumberOfPoints%, BaseSize%
  PRINT USING _
 "I found ### usable published distances in file \                        \"; _
                    I%, PUBDAT$ 
  PRINT "Please resolve this discrepency and start again!"
  CLOSE #2
  END
END IF

ON ERROR GOTO BadFile
DO
  NoGood$ = "F"
  PRINT:PRINT
  PRINT "What is the fully qualified name of the file"
  INPUT "   containing your Observations? : ", OBSDAT$
  IF OBSDAT$ = "" THEN
    PRINT:PRINT
    PRINT"I will take a null file name as a signal to stop executing!"
    CLOSE #2
    END
  END IF
  OPEN OBSDAT$ FOR INPUT AS #1
LOOP WHILE NoGood$ = "T"

ON ERROR GOTO 0

PRINT:PRINT
PRINT "Do you know the Manufacturer's Specifications"
INPUT "   for this instrument? (Y or N) : ", Answer$
WHILE Answer$ <> "N" AND Answer$ <> "n" AND Answer$ <> "Y" AND Answer$ <> "y"
  PRINT
  PRINT "I need a Y or N response.  Please try again!"
  PRINT:PRINT
  INPUT "Do you know the Manufacturer's Specifications? : ", Answer$
WEND

PPM = -999.999

IF Answer$ = "Y" OR Answer$ = "y" THEN
  PRINT:PRINT
  INPUT "What is the instrument's stated fixed error in mm? : " , Constant
  WHILE Constant < 0.0
    PRINT
    PRINT "I have a problem with that number.  Try again!"
    PRINT:PRINT
    INPUT "What is the instrument's stated fixed error in mm? : " , Constant
  WEND

  PRINT:PRINT
  INPUT "What is the instrument's stated scale error in PPM? : " , PPM
  WHILE PPM < 0.0
    PRINT
    PRINT "I have a problem with that number.  Try again!"
    PRINT:PRINT
    INPUT "What is the instrument's stated scale error in PPM? : " , PPM
  WEND

  IF PPM = 0.0 AND Constant = 0.0 THEN
    PRINT:PRINT
    PRINT "SORRY but you can not have a +-0/0 Manufacturer's specifications!"
    PRINT "Terminating execution!"
    CLOSE #1:CLOSE #2
    END
  END IF
END IF

PRINT:PRINT
PRINT "Please give me a line of Titling information"
PRINT "For example: WILD DI2000 (S/N 180010) ON CORBIN BASELINE"
PRINT
INPUT "TITLE is: ", Title$

CLS
PRINT #2,"      ";Title$
PRINT #2, USING "                     \        \"; DATE$
PRINT "      ";Title$
PRINT USING "                     \        \"; DATE$
PRINT:PRINT
PRINT #2," ":PRINT #2," "

IF PPM <> -999.999 THEN
  IF FIX(Constant) = Constant THEN
    PRINT USING "      Manufacturer's specifications: +### MM"; Constant
    PRINT #2, USING "      MANUFACTURER'S SPECIFICATIONS: +### MM"; Constant
  ELSE
    PRINT USING "      Manufacturer's specifications: +###.# MM"; Constant
    PRINT #2, USING "      MANUFACTURER'S SPECIFICATIONS: +###.# MM"; Constant
  END IF

  IF FIX(PPM) = PPM THEN
    PRINT USING "                                     +### PPM"; PPM
    PRINT #2, USING "                                     +### PPM"; PPM
  ELSE
    PRINT USING "                                     +###.# PPM"; PPM
    PRINT #2, USING "                                     +###.# PPM"; PPM
  END IF

  PRINT:PRINT
  PRINT #2," ":PRINT #2," "
END IF

ON ERROR GOTO ObsEOF1
WHILE NOT EOF(1)
  FromPoint = 0.0
  ToPoint = 0.0
  Observed = 0.0
  INPUT #1, FromPoint,ToPoint,Observed
  I% = 1
  DO
    IF (Published(I%,1) = FromPoint AND Published(I%,2) = ToPoint) _
        THEN EXIT LOOP
    IF (Published(I%,2) = FromPoint AND Published(I%,1) = ToPoint) _
        THEN EXIT LOOP
    I% = I% + 1
  LOOP WHILE I% <= BaseSize%

  IF I% > BaseSize% THEN
    PRINT
    PRINT "The following record from your User Observations File"
    PRINT "has an improper FROM or TO point number:"
    PRINT FromPoint "," ToPoint "," Observed
    PRINT #2,"The following record from your User Observations File"
    PRINT #2,"has an improper FROM or TO point number:"
    PRINT #2, FromPoint "," ToPoint "," Observed

  ELSEIF PPM <> -999.999 THEN
    Count = Count + 1
    DeltaMM = (Published(I%,3) - Observed) * 10^3
    Factor = Published(I%,3) / 1000.0
    Denominator = (Factor * PPM)^2 + Constant^2
    REG1 = REG1 + (Factor / Denominator)
    REG2 = REG2 + (Factor^2 / Denominator)
    REG3 = REG3 + (DeltaMM / Denominator)
    REG4 = REG4 + ((Factor * DeltaMM) / Denominator)
    REG5 = REG5 + (DeltaMM^2 / Denominator)
    REG6 = REG6 + (1.0 / Denominator)
  ELSE
    Count = Count + 1
    Delta(Count) =  Published(I%,3) - Observed
    SumPublished = SumPublished + Published(I%,3)
    SumPublishedSquared = SumPublishedSquared + Published(I%,3)^2
    SumDelta = SumDelta + Delta(Count)
    SumPublishedXDelta = SumPublishedXDelta + Published(I%,3) * Delta(Count)
  END IF
WEND
ON ERROR GOTO 0

IF Count < 3 THEN
  PRINT:PRINT
  PRINT USING _
  "SORRY \\there are only # good observations in \                        \";_
              "- " , Count , OBSDAT$
  PRINT "At least 3 observations are needed for this procedure"
  PRINT "Terminating execution!"
  CLOSE #1:CLOSE #2
  END

ELSEIF PPM <> -999.999 THEN
  REG0 = (REG6 * REG2) - (REG1^2)
  Scale = ( (REG6 * REG4) - (REG1 * REG3) ) / REG0
  REG8 = ( (REG2 * REG3) - (REG1 * REG4) ) / REG0
  Constant = REG8 / 1000.0
  VarUnitWeight = (REG5 - ( (REG4 * Scale) + (REG3 * REG8) ) ) / (Count - 2.0)
  StErrorScale = SQR( (REG6 / REG0) * VarUnitWeight )
  StErrorConstant = SQR( (REG2 / REG0) * VarUnitWeight)
  Covariance = (REG1 / REG0) * VarUnitWeight
  Correlation = Covariance / (StErrorScale * StErrorConstant)
  PRINT:PRINT #2," "
  PRINT:PRINT #2," "

  PRINT USING "    SCALE ERROR                   = +####.#### PPM";Scale
  PRINT USING "    STANDARD ERROR OF SCALE       = +####.#### PPM"; _
                   StErrorScale
  PRINT #2,USING "    SCALE ERROR                   = +####.#### PPM";Scale
  PRINT #2,USING "    STANDARD ERROR OF SCALE       = +####.#### PPM"; _
                      StErrorScale
  PRINT:PRINT #2," "

  PRINT USING "    CONSTANT ERROR                = +####.#### METERS";Constant
  PRINT USING "    STANDARD ERROR OF CONSTANT    = +####.#### METERS"; _
                   StErrorConstant / 1000.0
  PRINT #2,USING "    CONSTANT ERROR                = +####.#### METERS"; _
                      Constant
  PRINT #2,USING "    STANDARD ERROR OF CONSTANT    = +####.#### METERS"; _
                      StErrorConstant / 1000.0
  PRINT:PRINT #2, " "

  PRINT USING "    VARIANCE OF UNIT WEIGHT       = +####.####";VarUnitWeight
  PRINT #2,USING "    VARIANCE OF UNIT WEIGHT       = +####.####"; _
                      VarUnitWeight
  PRINT USING "    COVARIANCE                    = +####.####";Covariance
  PRINT #2,USING "    COVARIANCE                    = +####.####"; _
                      Covariance
  PRINT USING "    CORRELATION                   = +####.####";Correlation
  PRINT #2,USING "    CORRELATION                   = +####.####"; _
                      Correlation

ELSE
  MeanSumDelta = SumDelta / Count
  FOR I% = 1 to Count
    SumDeltaDiffSquared = SumDeltaDiffSquared + (Delta(I%) - MeanSumDelta)^2
  NEXT I%
  Scale = ( (Count * SumPublishedXDelta) - (SumPublished * SumDelta) ) / _
          ( (Count * SumPublishedSquared) - (SumPublished^2) )
  Constant = ( (SumPublishedSquared * SumDelta) - (SumPublished * _
              SumPublishedXDelta) ) / ( (Count * SumPublishedSquared) - _
              (SumPublished^2) )
  Sigma0Squared = ( SumDeltaDiffSquared - ( (Scale / Count) * ( (Count * _
                    SumPublishedXDelta) - (SumPublished * SumDelta) ) ) )_
                  / (Count - 2.0)
  StErrorScale = SQR( (Sigma0Squared * Count) / ( (Count * _
                    SumPublishedSquared) - (SumPublished^2) ) )
  ScaleStat = Scale / StErrorScale
  StErrorConstant = SQR ( (Sigma0Squared * SumPublishedSquared) / _
                        ( (Count * SumPublishedSquared) - (SumPublished^2) ) )
  ConstantStat = Constant / StErrorConstant
  PRINT:PRINT #2," "
  PRINT:PRINT #2," "

  PRINT USING "    SCALE ERROR                   = +###.#### PPM";Scale * 10^6
  PRINT USING "    STANDARD ERROR OF SCALE       = +###.#### PPM"; _
              StErrorScale * 10^6
  PRINT USING "    STUDENT'S t TEST              = +###.####";ScaleStat
  PRINT #2,USING "    SCALE ERROR                   = +###.#### PPM";Scale*10^6
  PRINT #2,USING "    STANDARD ERROR OF SCALE       = +###.#### PPM"; _
                 StErrorScale * 10^6
  PRINT #2,USING "    STUDENT'S t TEST              = +###.####";ScaleStat

  PRINT:PRINT #2," "

  PRINT USING "    CONSTANT ERROR                =   +#.#### METERS";Constant
  PRINT USING "    STANDARD ERROR OF CONSTANT    =   +#.#### METERS"; _
              StErrorConstant
  PRINT USING "    STUDENT'S t TEST              = +###.####";ConstantStat
  PRINT #2,USING "    CONSTANT ERROR                =   +#.#### METERS"; _
                 Constant
  PRINT #2,USING "    STANDARD ERROR OF CONSTANT    =   +#.#### METERS"; _
                 StErrorConstant
  PRINT #2,USING "    STUDENT'S t TEST              = +###.####";ConstantStat

END IF

CLOSE #1
OPEN OBSDAT$ FOR INPUT AS #1

IF PPM <> -999.999 THEN
  Scale = Scale / 10^6
END IF

PRINT #2," "
PRINT #2," "
PRINT #2," "
PRINT #2,"         PUBLISHED     OBSERVED   DIFFERENCE  RESIDUAL"
PRINT #2,"          (METERS)     (METERS)    (METERS)   (METERS)"
PRINT #2," "

ON ERROR GOTO ObsEOF2
WHILE NOT EOF(1)
  INPUT #1, FromPoint,ToPoint,Observed
  I% = 1
  DO
    IF (Published(I%,1) = FromPoint AND Published(I%,2) = ToPoint) _
        THEN EXIT LOOP
    IF (Published(I%,2) = FromPoint AND Published(I%,1) = ToPoint) _
        THEN EXIT LOOP
    I% = I% + 1
  LOOP WHILE I% <= BaseSize%

  IF I% <= BaseSize% THEN
    Difference = Published(I%,3) - Observed
    Residual = Difference - (Published(I%,3) * Scale) - Constant
    PRINT #2, USING "         ####.####    ####.####    +#.####    +#.####"; _
                     Published(I%,3),Observed,Difference,Residual
  END IF
WEND

PRINT:PRINT
PRINT "For a complete listing of the results, consult file"
PRINT "RESULTS.DAT in your presently connected directory."
CLOSE #1:CLOSE #2
END

BadFile:
  NoGood$ = "T"
  PRINT
  IF ERR = 75 OR ERR = 76 THEN
    PRINT "The path for this file is incorrect.  Try again!"
  ELSEIF ERR = 71 OR ERR = 72 THEN
    PRINT "The disk associated with this file is not ready.  Try again!"
  ELSE
    PRINT "The file was not found.  Try again!"
  END IF
  RESUME NEXT

PubEOF:
  IF ERR = 62 THEN
    PRINT
    PRINT "The last record in the Published Distances File does not have"
    PRINT "adequate data:"
    PRINT "   FROM Point = ", FIX(X1)
    PRINT "   TO Point   = ", FIX(X2)
    PRINT "   Distance   = ", X3
    PRINT
    PRINT "Please resolve this discrepency and start again!"
    CLOSE #1:CLOSE #2
    END
  ELSE
    RESUME NEXT
  END IF

ObsEOF1:
  IF ERR = 62 THEN
    PRINT
    PRINT "The last record in the User Observation File does not have"
    PRINT "adequate data:"
    PRINT "   FROM Point = ", FromPoint
    PRINT "   TO Point   = ", ToPoint
    PRINT "   Distance   = ", Observed
    PRINT #2," "
    PRINT #2,"The last record in the User Observation File does not have"
    PRINT #2,"adequate data:"
    PRINT #2,"   FROM Point = ", FromPoint
    PRINT #2,"   TO Point   = ", ToPoint
    PRINT #2,"   Distance   = ", Observed
  END IF
  RESUME NEXT

ObsEOF2:
  RESUME NEXT