MODULE PlanetGeometry;

IMPORT PGM, Math, LR3, LR3Extras;
IMPORT Wr, Rd, Text, Thread, Process, FileRd, OSError, ParseParams;
FROM Stdio IMPORT stdin, stdout, stderr;

PROCEDURE Main() =
  <* FATAL Wr.Failure, Thread.Alerted, Rd.Failure *>
  VAR totLight: LONG; 
      nHits: NAT;
      mkPGM: PGM.Reader;   (* Mask pixel reader *)
      mkMax: NAT;          (* Max pixel value in mask *)
      mkStep: LONG;        (* 1/mkMax *)
      mkC, otC: NAT;
  BEGIN
    WITH 
      d = GetData()
    DO
      CompleteFeatures(d);
      WriteData(stdout, d);
      Wr.Close(stdout);
    END
  END Main;
  
PROCEDURE FindCircle(
    READONLY x, y: ARRAY [0..2] OF LONG; 
    VAR ctrX, ctrY, radius: LONG;
  ) =
  (*
    Computes the center and radius of a circle, given 
    the coordinates of three points on its boundary.
  *)
  BEGIN
    WITH
      s0 = x[0]*x[0] + y[0]*y[0],
      s1 = x[1]*x[1] + y[1]*y[1],
      s2 = x[2]*x[2] + y[2]*y[2],
      
      A = + LR3Extras.Det(
        LR3.T{x[0], y[0], 1.0d0},
        LR3.T{x[1], y[1], 1.0d0},
        LR3.T{x[2], y[2], 1.0d0}
      ),
      
      D = - LR3Extras.Det(
        LR3.T{s0, y[0], 1.0d0},
        LR3.T{s1, y[1], 1.0d0},
        LR3.T{s2, y[2], 1.0d0}
      ),
      E = + LR3Extras.Det(
        LR3.T{s0, x[0], 1.0d0},
        LR3.T{s1, x[1], 1.0d0},
        LR3.T{s2, x[2], 1.0d0}
      ),
      F = - LR3Extras.Det(
        LR3.T{s0, x[0], y[0]},
        LR3.T{s1, x[1], y[1]},
        LR3.T{s2, x[2], y[2]}
      )
    DO
      ctrX := - 0.5d0 * D/A;
      ctrY := - 0.5d0 * E/A;
      radius := Math.sqrt(MAX(0.0d0, ctrX*ctrX + ctrY*ctrY - F/A))
    END
  END FindCircle;

PROCEDURE FindSunFromTerminator(
    READONLY x, y: ARRAY [0..1] OF LONG; 
    diskX, diskY, diskR: LONG;
    VAR sunDir: LR3.T;
  ) =
  (*
    Computes a sun direction vector, in screen coordinates, given the 
    disk center and radius, and two points "(x[i],y[i])" on
    the terminator.
  *)
  BEGIN
    WITH
      r2 = radius * radius,
      
      uX = x[0] - diskX,
      uY = y[0] - diskY,
      uZ = Math.sqrt(MAX(0.0d0, r2 - uX*uX - uY*uY)),
      
      vX = x[1] - diskX,
      vY = y[1] - diskY,
      vZ = Math.sqrt(MAX(0.0d0, r2 - vX*vX - vY*vY))
    DO
      (* Note that the screen coordinate system is left-handed... *)
      sunDir := LR3.Dir(LR3Extras.Cross(LR3.T{vX, vY, vZ}, LR3.T{uX, uY, uZ}))
    END
  END FindSun;
  
PROCEDURE FindMatrix(READONLY f: Features; VAR M: LR3x3) =
  VAR 
  BEGIN
    
  END FindMatrix;

PROCEDURE GetData (): Options =
  <* FATAL Thread.Alerted, Wr.Failure *>

  CONST
    MaxSize: NAT = 4096;
    MaxMaxVal: NAT = 256*256-1;
    MaxCoord: LONG = 10000.0d0;

  VAR d: Data;
  BEGIN
    WITH
      pp = NEW(ParseParams.T).init(stderr)
    DO

      TRY
        (* Get disk geometry: *)
        IF pp.keywordPresent("-limb") THEN
          VAR x, y: ARRAY [0..2] OF LONG;
          BEGIN
            FOR i := 0 TO 2 DO 
              x[i]:= pp.getNextLongReal(-MaxCoord, +MaxCoord);
              y[i]:= pp.getNextLongReal(-MaxCoord, +MaxCoord)
            END;
            FindCircle(x, y, d.diskX, d.diskY, d.diskR)
          END;
        ELSIF pp.keywordPresent("-disk") THEN
          d.diskX := pp.getNextLongReal(-MaxCoord, +MaxCoord);
          d.diskY := pp.getNextLongReal(-MaxCoord, +MaxCoord);
          d.diskR := pp.getNextLongReal(0.0d0, MaxCoord)
        ELSE
          pp.error("please specify either \"-limb\" or \"-disk\"")
        END;
        
        (* Get reference features: *)
        VAR NF: NAT := 0;
            f: REF Features := NEW(REF Features, 10);
        BEGIN
          WHILE pp.keywordPresent("-feature") DO
            IF NF >= NUMBER(f^) THEN
              WITH fNew = NEW(REF Features, 2*NF) DO
                SUBARRAY(fNew^, 0, NF) := f^;
                f := fNew
              END
            END;
            WITH fi = f[NF] DO
              fi.name := pp.getNext();
              IF pp.testNext("???") THEN
                <* ASSERT fi.testNext("???") *>
                fi.imgX := Huge; fi.imgY := Huge;
              ELSE
                fi.imgX := pp.getNextLongReal(-MaxCoord, +MaxCoord);
                fi.imgY := pp.getNextLongReal(-MaxCoord, +MaxCoord);
              END;
              IF pp.testNext("???") THEN
                <* ASSERT fi.testNext("???") *>
                fi.lat := Huge; fi.lng := Huge;
              ELSE
                fi.lat := pp.getNextLongReal(-90.0d0, +90.0d0);
                fi.lng := pp.getNextLongReal(-360.0d0, +360.0d0);
              END;
              IF pp.testNext("???") THEN
                <* ASSERT fi.testNext("???") *>
                <* ASSERT fi.testNext("???") *>
                fi.pos := LR3.Zero()
              ELSE
                fi.pos[0] := pp.getNextLongReal(-MaxCoord, +MaxCoord);
                fi.pos[1] := pp.getNextLongReal(-MaxCoord, +MaxCoord);
                fi.pos[2] := pp.getNextLongReal(-MaxCoord, +MaxCoord);
                fi.pos := LR3.Dir(fi.pos)
              END;
            END;
            INC(NF)
          END;
          d.f := NEW(REF Features, NF);
          d.f^ := SUBARRAY(f^, 0, NF)
        END;
          
        (* Get reference matrix: *)       
        IF pp.keywordPresent("-planetAxes") THEN
          d.M := GetMatrix(pp)
        ELSE
          FindMatrix(d.f, d.M)
        END;

        IF pp.keywordPresent("-sunDir") THEN
          d.sunDir := GetDir(pp)
        ELSIF pp.keywordPresent("-terminator") THEN
          VAR x, y: ARRAY [0..1] OF LONG;
          BEGIN
            FOR i := 0 TO 1 DO 
              x[i]:= pp.getNextLongReal(-0.5d0, FLOAT(d.NX, LONG)-0.5d0);
              y[i]:= pp.getNextLongReal(-0.5d0, FLOAT(d.NY, LONG)-0.5d0)
            END;
            FindSun(x, y, d.diskX, d.diskY, d.radius, d.sun)
          END;
        ELSE
          d.sun := LR3.T{0.0d0, 0.0d0, 1.0d0}
        END;
        
        IF pp.keywordPresent("-light") THEN
          WITH 
            maxLit = pp.getNextInt(1, d.maxVal)
          DO
            d.light := FLOAT(maxLit-1, LONG)/FLOAT(MAX(1, d.maxVal-1), LONG)
          END
        ELSE
          d.light := 1.0d0
        END;
        
        IF pp.keywordPresent("-dark") THEN
          WITH 
            maxLit = 1 + ROUND(d.light * FLOAT(d.maxVal-1, LONG)),
            minLit = pp.getNextInt(1, maxLit)
          DO
            d.dark := FLOAT(minLit-1, LONG)/FLOAT(MAX(1, d.maxVal-1), LONG)
          END
        ELSE
          d.dark := 1.0d0
        END;
        
        IF pp.keywordPresent("-bias") THEN
          d.bias := pp.getNextLongReal(-1.0d0, 1.0d0)
        ELSE
          d.bias := 0.0d0
        END;
        
        IF pp.keywordPresent("-mask") THEN
          d.mkRd := GetFileOption(pp);
        ELSE
          d.mkRd := NIL;
        END;
        

        pp.skipParsed();
        pp.finish();                                       
      EXCEPT                                                            
      | ParseParams.Error =>                                              
          Wr.PutText(stderr, "Usage: PlanetGeometry \\\n");
          Wr.PutText(stderr, "  -size <NX> <NY> [ -maxVal <M> ] \\\n");
          Wr.PutText(stderr, "  [ -mask <pgmfile> ] \\\n");
          Wr.PutText(stderr, "  [ -limb <X1> <Y1> <X2> <Y2> <X3> <Y3> | -disk <cX> <cY> <R> ] \\\n");
          Wr.PutText(stderr, "  [ -sun <X> <Y> <Z> | -terminator <X1> <Y1> <X2> <Y2> ] \\\n");
          Wr.PutText(stderr, "  [ -light <L> ] [-dark <K>] [ -bias <B> ]\n");
          Process.Exit(1);
      END;
    END;
    RETURN d
  END GetData;
  
PROCEDURE GetFileOption(pp: ParseParams.T): Rd.T RAISES {ParseParams.Error} =
  BEGIN
    WITH name = pp.getNext() DO
      IF Text.Equal(name, "-") THEN
        RETURN stdin
      ELSE
        TRY
          RETURN FileRd.Open(name)
        EXCEPT
        | OSError.E => 
            pp.error("can't open file \"" & name & "\"");
            RETURN NIL
        END
      END
    END
  END GetFileOption;

PROCEDURE Getmatrix(pp: ParseParams.T): LR3x3.T RAISES {} =
  VAR M: LR3x3.T;
  BEGIN
    FOR i := 0 TO 2 DO
      FOR j := 0 TO 2 DO
        M[i,j] := pp.getNextLongReal(-MaxCoord, +MaxCoord);
      END;
      M[i] := LR3.Dir(d.
    

  END~

          FOR i := 0 TO 2 DO d.sun[i] := pp.getNextLongReal(-1.0d0, 1.0d0) END;
          d.sun := LR3.Dir(d.sun);


BEGIN 
  Main();
END PlanetGeometry.