MODULE EvolutionDiagram EXPORTS Main;

(* 
  Writes an Encapsulated Postscript file showing the derivative
  field for a two-component reaction system. *)
  
IMPORT Wr, Thread, Process, PSPlot, ParseParams;
FROM Stdio IMPORT stderr;

TYPE 
  LONG = LONGREAL;

  Options = RECORD
      nSteps: CARDINAL;  (* Number of steps *)
    END;
    
PROCEDURE Main( ) =
  BEGIN
    WITH
      o = GetOptions()
    DO
      DrawEvolutionDiagram(
        name := "circ",
        DA := CircDA, DB := CircDB,
        vMin := -1.1d0, vMax := +1.1d0,
        dMax := +1.1d0,
        nSteps := o.nSteps
      );
    
      DrawEvolutionDiagram(
        name := "lotka",
        DA := LotkaDA, DB := LotkaDB,
        vMin := -1.0d0, vMax := +3.0d0,
        dMax := +2.5d0,
        nSteps := o.nSteps
      );
    END
  END Main;

TYPE
  Derivative = PROCEDURE (a, b: LONGREAL): LONGREAL;

PROCEDURE DrawEvolutionDiagram(
    name: TEXT; 
    DA, DB: Derivative;
    vMin, vMax: LONGREAL;
    dMax: LONGREAL;
    nSteps: CARDINAL;
  ) =
  <* FATAL Thread.Alerted, Wr.Failure *>

  PROCEDURE PutDir(f: PSPlot.File; x, y: INTEGER; dx, dy: LONGREAL)=
    BEGIN
      WITH
        x1 = FLOAT(x, LONGREAL) + 0.5d0,
        y1 = FLOAT(y, LONGREAL) + 0.5d0,
        x2 = x1 + dx,
        y2 = y1 + dy
      DO
        f.segment(x1, y1, x2, y2);
        f.circle(x1, y1, 0.3)
      END;
    END PutDir;

  BEGIN
    WITH
      nameEPS = name & "-ev-dir.eps",
      f = NEW(PSPlot.EPSFile).open(nameEPS),
      fn = FLOAT(nSteps, LONGREAL)
    DO
      Wr.PutText(stderr, "writing " & nameEPS & "...\n");
      f.setScale(PSPlot.Axis.X, 0.0d0, FLOAT(nSteps, LONG));
      f.setScale(PSPlot.Axis.Y, 0.0d0, FLOAT(nSteps, LONG));
      f.setFillColor(PSPlot.Black);
      FOR y := 0 TO nSteps - 1 DO 
        FOR x := 0 TO nSteps - 1 DO 
          WITH 
            a = vMin + (FLOAT(x, LONGREAL) + 0.5d0)/fn * (vMax - vMin),
            b = vMin + (FLOAT(y, LONGREAL) + 0.5d0)/fn * (vMax - vMin),
            da = DA(a, b)/dMax,
            db = DB(a, b)/dMax
          DO
            PutDir(f, x, y, da, db)
          END
        END
      END;
      f.frame();
      IF vMin< 0.0d0 AND vMax > 0.0d0 THEN 
        f.setLineWidth(0.1);
        f.setLineDashed(ARRAY OF REAL{2.0}, 1.0);
        f.coordLine(PSPlot.Axis.X, 0.0d0);
        f.coordLine(PSPlot.Axis.Y, 0.0d0);
      END;
      f.close();
    END;
  END DrawEvolutionDiagram;

PROCEDURE GetOptions (): Options =
  <* FATAL Thread.Alerted, Wr.Failure *>
  VAR o: Options;
  BEGIN
    WITH pp = NEW(ParseParams.T).init(stderr) DO                         
      TRY
        
        IF pp.keywordPresent("-nSteps") THEN
          o.nSteps := pp.getNextInt(0, 300);
        ELSE
          o.nSteps := 20
        END;

        pp.finish();                                       
      EXCEPT                                                            
      | ParseParams.Error =>                                              
          Wr.PutText(stderr, "Usage: EvolutionDiagram\\\n");
          Wr.PutText(stderr, "  [ -nSteps nn ] \n");        
          Process.Exit (1);                                              
      END;
    END;
    RETURN o
  END GetOptions;
  
PROCEDURE CircDA(a, b: LONGREAL): LONGREAL =
  BEGIN
    RETURN a * (1.0d0 - a*a - b*b)
  END CircDA;
  
PROCEDURE CircDB(a, b: LONGREAL): LONGREAL =
  BEGIN
    RETURN b * (1.0d0 - a*a - b*b)
  END CircDB;
  
PROCEDURE LotkaDA(a, b: LONGREAL): LONGREAL =
  BEGIN
    RETURN a * (1.0d0 - b)
  END LotkaDA;
  
PROCEDURE LotkaDB(a, b: LONGREAL): LONGREAL =
  BEGIN
    RETURN b * (a - 1.0d0)
  END LotkaDB;
  
BEGIN
  Main();
END EvolutionDiagram.