program Montef;

uses
   MyFuncs;
   
const
   CarP = 27.2753/365.2425;
   
   CH = 245; 
   CN = CH*7; { Fixed by Data, Number of Carrington rotations }
   CG = CH*2+1;
   
   BN = 18; { Fixed by Usoskin }
   CA = 7;
   CS = CN div CA;
   CR = CN*2;
   
   NCount = 100;
   
   CL = CR*NCount;
   SN = 100;


   NSile = 3; { Fixed by Data, approx. half a year }

   FlipFlop = True;
   Margins  = True; 
   Smoothing = True;
   Correlations = True;

   Smooth = 0.3;
   DeltaLow= 1.0/18;
   ProbDelta = 0.05;

   ChangeProb = 0.33333;
   StayProb = 0.33333;
   
   NextProb = (1.0-ChangeProb);
   DownProb = NextProb*(1.0-StayProb)/2.0;
   UpProb = 2*DownProb;
              
var
   Cntr,CLBase: Integer;
   Org1,Org2,Ph1,Ph2,Smooth1,Smooth2,Grid1,Grid2,
   SMax1,SMax2,Scat,ScatSum,Slope: array of Double;


function NewPhase(Prob: Double): Double;

begin
   Result:=Random;
   if Correlations and (Prob>=0.0) then begin
      if Result<DownProb then
         Result:=Prob-ProbDelta
      else if Result<UpProb then   
         Result:=Prob+ProbDelta
      else if Result<NextProb then
         Result:=Prob
      else
         Result:=Random;   
      if Result<0.0 then
         Result:=Result+1.0
      else if Result>1.0 then
         Result:=Result-1.0;  
   end;
end;

function PhaseDiff(p1,p2: Double): Boolean;

begin
   if not Smoothing then
      Result:=True
   else begin   
      if p1>=p2 then begin
         if ((p1-p2)<=Smooth) or
            ((p2+1.0-p1)<=Smooth) then
            Result:=False
         else
            Result:=True;   
      end else begin
         if ((p2-p1)<=Smooth) or
            ((p1+1.0-p2)<=Smooth) then
            Result:=False
         else
            Result:=True;         
      end;   
   end;
end;

procedure Generate;

var
   CurProb1,CurProb2: Double;
   i: Integer;
   OK: Boolean;
   
begin
   CurProb1:=-1.0;
   CurProb2:=-1.0;
   for i:=0 to CN-1 do begin
      OK:=False;
      repeat
         CurProb1:=NewPhase(CurProb1);
         CurProb2:=NewPhase(CurProb2);
         if PhaseDiff(CurProb1,CurProb2) then begin
            OK:=True;
            Ph1[i]:=CurProb1;
            Ph2[i]:=CurProb2;
         end;
      until OK;   
   end;
   for i:=0 to CN-1 do begin
      Org1[i]:=Ph1[i];
      Org2[i]:=Ph2[i];
   end;
end;

procedure CheckCrossings;

var
   DeltaHigh,Cur1,Cur2,Nxt1,Nxt2,gmin: Double;
   Cross1,Cross2: Boolean;
   i,j,jmin: Integer;
   Gains: array [0..7] of Double;

function Gain(x1,x2: Double): Double;

var
   Def1,Def2: Double;
   
begin
   Def1:=Cur1-Frac(Cur1)+x1;
   Def2:=Cur2-Frac(Cur2)+x2;
   if Abs(Def1-Def2)<1.0 then
      Result:=Abs(Def1-Cur1)+Abs(Def2-Cur2)
   else
      Result:=100;   
end;

procedure Put(x1,x2: Double);

var
   Def1,Def2: Double;
   
begin
   Def1:=Cur1-Frac(Cur1)+x1;
   Def2:=Cur2-Frac(Cur2)+x2;
   Ph1[i+1]:=Def1;
   Ph2[i+1]:=Def2;   
end;
   
begin
   if (DeltaLow<=0) or (DeltaLow>=1.0) then
      Exit; 

   DeltaHigh:=1.0-DeltaLow;

   for i:=0 to CN-2 do begin
      Cur1:=ph1[i];
      Cur2:=ph2[i];
      Nxt1:=ph1[i+1];
      Nxt2:=ph2[i+1];

      Cross1:=(Cur1>=DeltaHigh) and (Nxt1<=DeltaLow);
      Cross2:=(Cur2>=DeltaHigh) and (Nxt2<=DeltaLow);
         
      for j:=0 to 7 do
         Gains[j]:=100;
               
      Gains[0]:=Gain(Nxt1,Nxt2);
      Gains[1]:=Gain(Nxt2,Nxt1);
         
      if Cross1 then begin
         Gains[2]:=Gain(Nxt1+1.0,Nxt2);
         Gains[3]:=Gain(Nxt2,Nxt1+1.0);
      end;
         
      if Cross2 then begin
         Gains[4]:=Gain(Nxt1,Nxt2+1.0);
         Gains[5]:=Gain(Nxt2+1.0,Nxt1);
      end;
         
{      if Cross2 and Cross1 then begin
         Gains[6]:=Gain(Nxt1+1.0,Nxt2+1.0);
         Gains[7]:=Gain(Nxt2+1.0,Nxt1+1.0);
      end;
}         
      gmin:=Gains[0];
      jmin:=0;

      for j:=1 to 7 do begin
         if Gains[j]<gmin then begin
            gmin:=Gains[j];
            jmin:=j;
         end;    
      end;      

      case jmin of
      0: Put(Nxt1,Nxt2);
      1: Put(Nxt2,Nxt1);
      2: Put(Nxt1+1.0,Nxt2);
      3: Put(Nxt2,Nxt1+1.0);
      4: Put(Nxt1,Nxt2+1.0);
      5: Put(Nxt2+1.0,Nxt1);
      6: {Put(Nxt1+1.0,Nxt2+1.0)}; 
      7: {Put(Nxt2+1.0,Nxt1+1.0)};
      end;
   end;
end;

procedure DoFlipFlop;

var
   Tmp: Double;
   i: Integer;
   
begin
   for i:=0 to CN-1 do
      if Ph1[i]>Ph2[i] then begin
         Tmp:=Ph1[i];
         Ph1[i]:=Ph2[i];
         Ph2[i]:=Tmp;
      end;
end;

procedure ReportPhases;

var
   i: Integer;
   
begin
   WriteLn('c Hit1 Hit2');
   for i:=0 to CN-1 do
      WriteLn(i*CarP:20:6,' ',Ph1[i]:12:6,' ', Ph2[i]:12:6);
end;

procedure SmoothPhases;

var
   MSile,j,i,k,iarg: Integer;
   Sum1,Sum2,Arg,Fg: Double;
   
begin
   MSile:=NSile*2+1;
   i:=0;
   for j:=0 to CH-1 do begin
      Sum1:=0.0;
      Sum2:=0.0;
      for k:=1 to MSile do begin
         Sum1:=Sum1+Ph1[i];
         Sum2:=Sum2+Ph2[i];
         Inc(i);
      end;
      Smooth1[j]:=Sum1/MSile;
      Smooth2[j]:=Sum2/MSile;
   end;
   j:=1;
   for i:=0 to CH-2 do begin
      Grid1[j]:=Smooth1[i];
      Grid2[j]:=Smooth2[i];
      Inc(j);
      Grid1[j]:=0.5*(Smooth1[i]+Smooth1[i+1]);
      Grid2[j]:=0.5*(Smooth2[i]+Smooth2[i+1]);
      Inc(j);
   end;
   Grid1[0]:=Smooth1[0]-0.5*(Smooth1[1]-Smooth1[0]);
   Grid2[0]:=Smooth2[0]-0.5*(Smooth2[1]-Smooth2[0]);
   Grid1[CG-2]:=Smooth1[CH-1];
   Grid2[CG-2]:=Smooth2[CH-1];
   Grid1[CG-1]:=Smooth1[CH-1]-0.5*(Smooth1[CH-2]-Smooth1[CH-1]);
   Grid2[CG-1]:=Smooth2[CH-1]-0.5*(Smooth2[CH-2]-Smooth2[CH-1]);
   for i:=0 to CN-1 do begin
      Arg:=(CG-1)*(i+0.5)/CN;
      Fg:=Frac(Arg);
      iarg:=Trunc(Arg);
      SMax1[i]:=(1.0-Fg)*Grid1[iarg]+Fg*Grid1[iarg+1];
      SMax2[i]:=(1.0-Fg)*Grid2[iarg]+Fg*Grid2[iarg+1];
   end;   
end;

procedure ReportGrids;

var
   i: Integer;
   
begin
   WriteLn('c Grid1 Grid2');
   for i:=0 to CG-1 do
      WriteLn(i,' ',Grid1[i]:20:6,' ',Grid2[i]:20:6);
end;

procedure ReportTrends;

var
   i: Integer;
   
begin
   WriteLn('c Org1 Org2');
   for i:=0 to CN-1 do begin
      if Abs(Frac(Ph1[i])-Org1[i])<1.0e-6 then
         WriteLn(i*CarP:20:6,' ',Ph1[i],' ',Ph2[i])
      else
         WriteLn(i*CarP:20:6,' ',Ph2[i],' ',Ph1[i]);   
   end;
end;

procedure ReportMeans;

var
   i: Integer;

begin
   WriteLn('c SMax1 SMax2');
   for i:=0 to CN-1 do
      WriteLn(i*CarP:20:6,' ',SMax1[i]:12:6,' ',SMax2[i]:12:6);
end;

procedure ComputeScatter;

var
   i: Integer;
   
begin
   for i:=0 to CN-1 do begin
      Scat[i*2]:=Ph1[i]-SMax1[i];
      Scat[i*2+1]:=Ph2[i]-SMax1[i];
   end;
end;

procedure ReportScatter;

var
   i: Integer;
   
begin
   WriteLn('c Sc');
   for i:=0 to CN-1 do begin
      WriteLn(2*i*CarP:20:6,' ',Scat[2*i]:12:6);
      WriteLn((2*i+1)*CarP:20:6,' ',Scat[2*i+1]:12:6);
   end;
   WriteLn('c Sing1 Sing2');
   for i:=0 to CN-1 do
      WriteLn(i*CarP:20:6,' ',Scat[2*i]:12:6,' ',Scat[2*i+1]:12:6);
end;

procedure AccumulateScatter;

var
   i: Integer;
   
begin
   for i:=0 to CR-1 do
      ScatSum[CLBase+i]:=Scat[i];
   CLBase:=CLBase+CR;     
end;

procedure RegisterSlope;

begin
   Slope[Cntr]:=0.5*Abs(Ph1[0]+Ph2[0]-Ph1[CN-1]-Ph2[CN-1]);
end;

procedure ReportSummary;

var
   i: Integer;
   
begin
   WriteLn('c Slope');
   for i:=0 to NCount-1 do
      WriteLn(i,' ',Slope[i]); 
   WriteLn('c Scat');
   for i:=0 to CL-1 do
      WriteLn(i,' ',ScatSum[i]);
end;

      
begin
   SetLength(Ph1,CN);
   SetLength(Ph2,CN);
   SetLength(Org1,CN);
   SetLength(Org2,CN);
   SetLength(Smooth1,CH);
   SetLength(Smooth2,CH);
   SetLength(Grid1,CG);
   SetLength(Grid2,CG);
   SetLength(SMax1,CN);
   SetLength(SMax2,CN);
   SetLength(Scat,CR);
   Randomize;
   Generate;
   if Margins then   
      CheckCrossings;
   if FlipFlop then
      DoFlipFlop;
   ReportPhases;
   SmoothPhases;
   ReportGrids;
   ReportTrends;
   ReportMeans;
   ComputeScatter;
   ReportScatter;
   if NCount>1 then begin
      SetLength(ScatSum,CL);
      SetLength(Slope,NCount);
      CLBase:=0;
      for Cntr:=0 to NCount-1 do begin
         Generate;
         if Margins then
            CheckCrossings;
         if FlipFlop then
            DoFlipFlop;
         RegisterSlope;   
         SmoothPhases;
         ComputeScatter;
         AccumulateScatter;
       end;
       ReportSummary;
   end;   
end.