Unit Utils;

interface

uses
   SysUtils,
   Graphics;

const
   Max_Rotations = 2500;
   Max_Data = 20000;
   Max_Source = 300000;
   Max_Records = 70;

const
   CR_Period = 27.2753;
   T_C = 25.38;

type
   TMethod = (mtHits,mtSimple,mtDouble,mtShifts);
   TGridEvaluator = function: Double;
   TGridEvaluatorProfile=function(Deg: Double): Double;
   TGrid = array of array of array of Single;
   TSlice = array of array of Single;
   TVector = array of Single;

var
   Evaluator: TGridEvaluator;
   EvaluatorProfile: TGridEvaluatorProfile;

   rmin,rmax,ldim,odim,bdim: Integer;

   CR: array [1..Max_Rotations] of Boolean;
   Mean: array [1..Max_Rotations] of Double;
   Correction: array [1..Max_Rotations] of Double;
   Shift: array [1..Max_Rotations] of Double;

   Cnt: array [1..2500] of Integer;
   Area: array [1..Max_Rotations,1..Max_Records] of Single;
   Longitude: array [1..Max_Rotations,1..Max_Records] of Single;

   Criterion: TGrid;
   Slice: TSlice;
   Vector: TVector;
   
   Min_Rot,Max_Rot: Integer;
   Min_Lat,Max_Lat: Double;
   JD_Min,JD_Max: Double;
   Min_Year,Max_Year: Integer;
   Min_Area,Max_Area: Integer;
   Is_Corrected: Boolean;
   Include_Singles: Boolean;


   Rot_Count: Integer;
   Hit_D: Integer;
   Norm_Coef: Double;

procedure Compile_Full_Dataset(MinR,MaxR: Integer; IsNorth: Boolean; Name: String);
procedure Compile_Strip(MinR,MaxR: Integer; IsNorth: Boolean; LMin,LMax: Double; Name: String);
procedure Compute_Method(LambdaMin,LambdaMax,LambdaStep,
                         OmegaMin,OmegaMax,OmegaStep,
                         BMin,BMax,BStep: Double; 
                         Name: String; 
                         Method: TMethod; PhysMode,Repar: Boolean; 
                         var Lambda,Omega,B: Double);
procedure Write_Shifts(Lambda,Omega,B: Double; Name,ResName: String; PhysMode: Boolean);
procedure Write_Corrections(Name,ResName: String);
procedure Write_Results(Lambda,Omega,B: Double; Name,ResName: String);

procedure RestoreDefaults;

procedure ReadDataFromSource(Source,Data: String; Is_North: Boolean);

procedure ReadMeansFromData(Data: String);
procedure InterpolateMeans;
procedure WriteMeans(Curve: String);

procedure ReadMeans(Curve: String);
procedure SmoothMeans(M: Integer);

procedure ComputeCorrections;
function SlopeOfCorrections: Double;
function ReparametrizeCorrections: Double;
procedure WriteCorrections(Crs: String);

procedure ComputeShifts(Lambda,Omega,B: Double; PhysMode: Boolean);
procedure WriteShifts(Shs: String);

procedure ReadData(Data: String; Use_Longitudes: Boolean);
procedure CheckData(Renorm_Areas: Boolean);
function CountRotations: Integer;
function CountRecords: Integer;
function ComputeNorm: Double;



procedure Compute_Grid(LambdaMin,LambdaMax,LambdaStep,OmegaMin,OmegaMax,OmegaStep,BMin,BMax,BStep: Double;PhysMode: Boolean);
procedure Compute_Omega_B_Slice;
procedure Compute_Omega_Vector;
procedure Compute_B_Vector;
procedure Compute_Lambda_Vector;
function Compute_Grid_Minimum(var Lambda,Omega,B: Double): Double;

procedure Paint_Slice(Name: String);
procedure Write_Vector(Name: String);

procedure Write_Profile(Name: String);
procedure Write_Phases(Name: String);
procedure Write_Distribution(Name: String);

procedure Eliminate_Longitudes;

function Compute_Hits: Double;
function Compute_Simple: Double;
function Compute_From_Shifts: Double;
function Compute_NAxi: Double;
function Compute_Double: Double;

implementation

var
   Grid_LambdaMin,
   Grid_LambdaStep,
   Grid_OmegaMin,
   Grid_OmegaStep,
   Grid_BMin,
   Grid_BStep,
   Vector_Min,
   Vector_Step: Double;

function ReadInt(S: String; st,cnt: Integer): Integer;

var
   D: String;
   V,Code: Integer;

begin
   D:=Trim(Copy(S,st,cnt));
   Val(D,V,Code);
   if Code=0 then
      Result:=V
   else
      Result:=0;   
end;

function ReadFloat(S: String; st,cnt: Integer): Double;

var
   D: String;
   V: Double;   
   Code: Integer;

begin
   D:=Trim(Copy(S,st,cnt));
   Val(D,V,Code);
   if Code=0 then
      Result:=V
   else
      Result:=0.0;   
end;

function JD(Day,Month,Year:Integer;Part:Double):Double;

var 
   A: Double;
   B: Integer;

begin
   A:=10000.0*Year+100.0*Month+Day;
   if (Month<=2) then begin 
      Month:=Month+12; 
      Year:=Year-1;
   end;
   if (A<=15821004.1) then
      B:=-2+trunc((Year+4716)/4)-1179
   else 
      B:=trunc(Year/400)-trunc(Year/100)+trunc(Year/4);
   A:=365.0*Year-679004.0;
   Result:=2400000.5+A+B+trunc(30.6001*(Month+1))+Day+Part;
end;

procedure ReadDataFromSource(Source,Data: String; Is_North: Boolean);

var
   F,G: Text;
   JD0: Double;
   Hit,Hit2: array of Boolean;
   i,Year,Month,Day,Cri,Nbr,Def,OldCri,isum,Cd,Area: Integer; 
   Part,Julius,CR,Longitude,Latitude,Meridian,Sum,Norm,Mean: Double;
   DoIt: Boolean;
   S: String;
 
begin
   Assign(F,Source);
   Reset(F);
   Assign(G,Data);
   Rewrite(G);
   JD0:=JD(9,11,1853,0.000);
   WriteLn(G,'r Zero Epoch = ',JD0:12:4);
   SetLength(Hit,Max_Source);
   SetLength(Hit2,Max_Data);
   for i:=0 to Max_Source-1 do
      Hit[i]:=True;
   for i:=0 to Max_Data-1 do
      Hit2[i]:=True;

   WriteLn(G,'c CR Y M D P N Def Area Long Lat Mean');
   WriteLn(G,'l -100000.0');
   WriteLn(G,'h 100000.0'); 

   OldCri:=0;

   Sum:=0.0;
   Norm:=0.0;
   isum:=0;  

   while not eof(F) do begin
      ReadLn(F,S);
      Year:=ReadInt(S,1,4);
      Month:=ReadInt(S,5,2);
      Day:=ReadInt(S,7,2);
      Part:=ReadFloat(S,9,4);
      Julius:=JD(Day,Month,Year,Part);
      CR:=Trunc((Julius-JD0)/CR_Period)+1;
      Cri:=Round(CR);
      Nbr:=ReadInt(S,13,8);

      if Is_Corrected then
         Area:=ReadInt(S,40,5)
      else
         Area:=ReadInt(S,30,5);

      Longitude:=ReadFloat(S,58,5);
      Latitude:=ReadFloat(S,64,5);
      Meridian:=ReadFloat(S,70,5);

      DoIt:=False;
      Def:=0;
          
      if Year>=1977 then begin
         if Hit2[Nbr] then
            DoIt:=True;
         Hit2[Nbr]:=False;
      end else begin
         if Nbr<=23738 then begin

            { Group entry }
            { Some numbers are repeated, check for 21-th column }
 
            Def:=ReadInt(S,21,1);
            Cd:=Nbr*10+Def;
            if Hit[Cd] then
               DoIt:=True;
            Hit[Cd]:=False;
         end else if Include_Singles then
            DoIt:=True;
      end;

      { Sanity Check }

      if Longitude>360.0 then begin
         if Longitude<370.0 then
            Longitude:=Longitude-360.0
         else
            DoIt:=False;
      end; 

      { Filters }

      if (Year<Min_Year) or (Year>Max_Year) then
         DoIt:=False;

      if Is_North then begin
         if Latitude<0.0 then
            DoIt:=False; 
      end else begin
         if Latitude>0.0 then
            DoIt:=False
         else
            Latitude:=-Latitude;
      end; 

      if (Area<Min_Area) or (Area>Max_Area) then
         DoIt:=False;

      if (Latitude<Min_Lat) or (Latitude>Max_Lat) then
         DoIt:=False; 

      if (Julius<JD_Min) or (Julius>JD_Max) then
         DoIt:=False; 

      if (Cri<Min_Rot) or (Cri>Max_Rot) then 
         DoIt:=False; 

      if Meridian<-1000.0 then
         DoIt:=False;

      if DoIt then begin
         if Cri<>OldCri then begin
             if (isum=0) or (Abs(Norm)<1.0e-10) then
                Mean:=0.0
             else
                Mean:=Sum/Norm;
             Sum:=0.0;
             Norm:=0.0;
             isum:=0;
         end else
             Mean:=0.0;
         if OldCri<>0 then 
            WriteLn(G,Mean:12:8);
         Sum:=Sum+Area*Latitude;
         Norm:=Norm+Area;
         isum:=isum+1;
         Write(G,Julius:9:4,' ',CR:4:0,' ',Year:4,' ',Month:2,' ',Day:2,' ',Part:6:4,' ',
                 Nbr:8,' ',Def:2,' ',Area:6,' ',Longitude:6:1,' ',Latitude:6:1,' ');
         OldCri:=Cri;
      end;
   end;

   if (isum=0) or (Abs(Norm)<1.0e-10) then
      Mean:=0.0
   else begin
      Mean:=Sum/Norm;
   end;

   WriteLn(G,Mean:12:8);  
   Close(F);
   Close(G);
end;

procedure ReadMeansFromData(Data: String);

var
   i: Integer;
   F: Text;
   S: String;

begin
   for i:=1 to Max_Rotations do begin
      CR[i]:=False;
      Mean[i]:=0.0;
   end;

   rmin:=Max_Rotations+1;
   rmax:=0;

   Assign(F,Data);
   Reset(F);
   while not eof(F) do begin
      ReadLn(F,S);
      if S[2]=' ' then
         Continue;
      i:=ReadInt(S,14,4);
      if i<>0 then begin
         Mean[i]:=ReadFloat(S,70,8);
         CR[i]:=True;
         if i<rmin then
            rmin:=i;
         if i>rmax then
            rmax:=i; 
      end;
   end;
   Close(F);
end;

procedure InterpolateMeans;

var
   i,j,i1,i2: Integer;
   left,right,Step,Ipv: Double;

begin
   for i:=rmin to rmax do
      if (not CR[i]) or (Mean[i]=0.0) then begin
         i1:=i;
         i2:=i;
         while i2<rmax do begin
            if CR[i2+1] and (Mean[i2+1]<>0.0) then
               Break;
            Inc(i2);   
         end;
         if i1=rmin then begin
            right:=Mean[i2+1];
            left:=right; 
         end else if i2=rmax then begin
            left:=Mean[i1-1];
            right:=Left;  
         end else begin
            left:=Mean[i1-1];
            right:=Mean[i2+1]; 
         end;
         Step:=(right-left)/(i2-i1+2);
         Ipv:=Left+Step;
         for j:=i1 to i2 do begin
             CR[j]:=True;
             Mean[j]:=Ipv;
             Ipv:=Ipv+Step; 
         end;
      end;
end;

procedure WriteMeans(Curve: String);

var
   i: Integer;
   F: Text;

begin
   Assign(F,Curve+'.dat');
   Rewrite(F);
   WriteLn(F,'c Mean');
   for i:=rmin to rmax do
      if CR[i] then
         WriteLn(F,i:4,' ',Mean[i]:14:8);
   Close(F);  
end;
  
procedure ReadMeans(Curve: String);

var
   i: Integer;
   F: Text;
   S: String;

begin
   for i:=1 to Max_Rotations do begin
      CR[i]:=False;
      Mean[i]:=0.0;
   end;

   rmin:=Max_Rotations+1;
   rmax:=0;

   Assign(F,Curve);
   Reset(F);
   while not eof(F) do begin
      ReadLn(F,S);
      if S[2]=' ' then
         Continue;
      i:=ReadInt(S,1,4);
      if i<>0 then begin
         Mean[i]:=ReadFloat(S,6,8);
         CR[i]:=True;
         if i<rmin then
            rmin:=i;
         if i>rmax then
            rmax:=i; 
      end;
   end;
   Close(F);
end;

procedure SmoothMeans(M: Integer);

var
   i,j,j1,j2,k: Integer;
   Cn: Double;
   Back: array [1..Max_Rotations] of Double;
 
begin
   for i:=1 to Max_Rotations do
      if CR[i] then
         Back[i]:=Mean[i]
      else
         Back[i]:=0.0;
        
   for i:=rmin to rmax do 
      if CR[i] then begin
         j1:=i-M;
         j2:=i+M;
         k:=0;
         Cn:=0.0;
         for j:=j1 to j2 do
            if (j>1) and (j<Max_Rotations) then begin
               Inc(k);
               Cn:=Cn+Back[j] 
            end;   
         Mean[i]:=Cn/k;
      end;
end;

procedure ComputeCorrections;

var
   i: Integer;
   Cum: Double;

begin
   for i:=rmin to rmax do
      Correction[i]:=Sqr(Sin(2*Pi*Mean[i]/360));
   Cum:=0.0;
   for i:=rmin to rmax do begin
      Cum:=Cum+Correction[i];
      Correction[i]:=Cum;
   end;
end;

function SlopeOfCorrections: Double;

begin
   Result:=(Correction[rmax]-Correction[rmin])/(rmax-rmin);
end;

function ReparametrizeCorrections: Double;

var
   i,idef: Integer;

begin
   idef:=rmin-1;
   Result:=SlopeOfCorrections;
   for i:=rmin to rmax do
      Correction[i]:=Correction[i]-Result*(i-idef);
end;


procedure WriteCorrections(Crs: String);

var
   F: Text;
   i: Integer;
   Cn: Double;

begin
   Assign(F,Crs+'.dat');
   Rewrite(F);
   WriteLn(F,'c Corrs');
   WriteLn(F,'l -100000');
   WriteLn(F,'h  100000');
   for i:=rmin to rmax do begin
      Cn:=Correction[i]; 
      WriteLn(F,i:4,' ',Cn:10:6);
   end;
   Close(F); 
end;
 
procedure ComputeShifts(Lambda,Omega,B: Double; PhysMode: Boolean);

var
   i,idef: Integer;
   Omega_C,Df: Double;

begin
   idef:=rmin-1;
   Omega_C:=360.0/T_C;
   if PhysMode then begin
      Df:=Omega-Omega_C;
      for i:=rmin to rmax do
         Shift[i]:=Lambda+T_C*(Df*(i-idef)-B*Correction[i]);
   end else begin
      Df:=Omega_C-Omega;
      for i:=rmin to rmax do
         Shift[i]:=Lambda+T_C*(Df*(i-idef)+B*Correction[i]);
   end;
end;

procedure WriteShifts(Shs: String);

var
   F: Text;
   i: Integer;
   Cn: Double;

begin
   Assign(F,Shs+'.dat');
   Rewrite(F);
   WriteLn(F,'l -100000');
   WriteLn(F,'h  100000');
   WriteLn(F,'c Shifts');
   for i:=rmin to rmax do 
      if CR[i] then begin
         Cn:=Shift[i]; 
         WriteLn(F,i:4,' ',Cn:10:6);
      end;
   Close(F); 
end;

procedure ReadData(Data: String; Use_Longitudes: Boolean);

var
   F: Text;
   i: Integer;
   S: String;
   Lo,Ar: Double;

begin
   for i:=1 to Max_Rotations do
      Cnt[i]:=0;
   Assign(F,Data);
   Reset(F);
   while not eof(F) do begin
      ReadLn(F,S);
      if S[2]=' ' then
          Continue;
      i:=ReadInt(S,14,4);
      if CR[i] then begin
          Lo:=ReadFloat(S,56,6);
          Ar:=ReadInt(S,49,6);
          Cnt[i]:=Cnt[i]+1;
          if Cnt[i]<=Max_Records then begin
             Area[i,Cnt[i]]:=Ar;
             if Use_Longitudes then
                Longitude[i,Cnt[i]]:=Lo
             else
                Longitude[i,Cnt[i]]:=0.0;
          end;
      end; 
   end;
   Close(F);
end;

procedure CheckData(Renorm_Areas: Boolean);

var
   i,j,k: Integer;
   Sum: Double;

begin
   for i:=rmin to rmax do
      if CR[i] then begin
         if Cnt[i]=0 then
            CR[i]:=False
         else begin
            Sum:=0.0;
            for j:=1 to Cnt[i] do
               Sum:=Sum+Area[i,j];
            if Sum=0.0 then
               CR[i]:=False
            else begin

               { Ignore spots with zero area }

               k:=0;
               for j:=1 to Cnt[i] do 
                  if Area[i,j]>0.0 then begin
                     Inc(k);
                     Area[i,k]:=Area[i,j];
                     Longitude[i,k]:=Longitude[i,j];
                  end; 
               Cnt[i]:=k;
               if Renorm_Areas then
                  for j:=1 to Cnt[i] do 
                     Area[i,j]:=Area[i,j]/Sum; 
            end;
         end;
      end;
end;

function CountRotations: Integer;

var
   i: Integer;

begin
   Result:=0;
   for i:=rmin to rmax do
      if CR[i] then
         Inc(Result); 
end;

function CountRecords: Integer;

var
   i: Integer;

begin
   Result:=0;
   for i:=rmin to rmax do
      if CR[i] then
         Result:=Result+Cnt[i]
end;

function ComputeNorm: Double;

var
   i,j: Integer;

begin
   Result:=0.0;
   for i:=rmin to rmax do
      if CR[i] then 
         for j:=1 to Cnt[i] do
            Result:=Result+Area[i,j];
end;

procedure Compute_Grid(LambdaMin,LambdaMax,LambdaStep,OmegaMin,OmegaMax,OmegaStep,BMin,BMax,BStep: Double; PhysMode: Boolean);

var
   i,j,k: Integer;

begin
   if (LambdaStep=0.0) or (LambdaMax<LambdaMin) then
      ldim:=1
   else begin
      ldim:=Round((LambdaMax-LambdaMin)/LambdaStep);
      if (LambdaMin+(ldim-1)*LambdaStep)<LambdaMax then
         Inc(ldim); 
   end;
   if (OmegaStep=0.0) or (OmegaMax<OmegaMin) then
      odim:=1
   else begin
      odim:=Round((OmegaMax-OmegaMin)/OmegaStep);
      if (OmegaMin+(odim-1)*OmegaStep)<OmegaMax then
         Inc(odim);
   end; 
   if (BStep=0.0) or (BMax<BMin) then
      bdim:=1
   else begin
      bdim:=Round((BMax-BMin)/BStep);
      if (BMin+(odim-1)*BStep)<BMax then
         Inc(bdim);
   end; 
   Grid_LambdaMin:=LambdaMin;
   Grid_LambdaStep:=LambdaStep;
   Grid_OmegaMin:=OmegaMin;
   Grid_OmegaStep:=OmegaStep;
   Grid_BMin:=BMin;
   Grid_BStep:=BStep;
   SetLength(Criterion,ldim,odim,bdim);
   for i:=0 to ldim-1 do
      for j:=0 to odim-1 do
          for k:=0 to bdim-1 do begin
              ComputeShifts(LambdaMin+i*LambdaStep,OmegaMin+j*OmegaStep,BMin+k*BStep,PhysMode);
              Criterion[i,j,k]:=Evaluator;
          end;
   WriteLn('r ');
   WriteLn('r Lambda Min   = ',LambdaMin:12:4);
   WriteLn('r Lambda Max   = ',LambdaMax:12:4);
   WriteLn('r Lambda Step  = ',LambdaStep:12:4);
   WriteLn('r Lambda Count = ',ldim:6); 
   WriteLn('r Omega  Min   = ',OmegaMin:12:8);
   WriteLn('r Omega  Max   = ',OmegaMax:12:8);
   WriteLn('r Omega  Step  = ',OmegaStep:12:8);
   WriteLn('r Omega  Count = ',odim:6);
   WriteLn('r B      Min   = ',BMin:12:6);
   WriteLn('r B      Max   = ',BMax:12:6);
   WriteLn('r B      Step  = ',BStep:12:6);
   WriteLn('r B      Count = ',bdim:6);
   WriteLn('r Total  Count = ',(ldim*odim*bdim):8);  
   WriteLn('r '); 
end;

procedure Compute_Omega_B_Slice;

var
   i,j,k: Integer;
   cmin: Double;

begin
   SetLength(Slice,odim,bdim);
   for j:=0 to odim-1 do
      for k:=0 to bdim-1 do begin
         cmin:=Criterion[0,j,k];
         for i:=1 to ldim-1 do
            if Criterion[i,j,k]<cmin then 
               cmin:=Criterion[i,j,k];
         Slice[j,k]:=cmin;
      end;    
end;

function Compute_Grid_Minimum(var Lambda,Omega,B: Double): Double;

var
   i,j,k,imin,jmin,kmin: Integer;
   cmin: Double;

begin
   cmin:=Criterion[0,0,0];
   imin:=0;
   jmin:=0;
   kmin:=0;
   for i:=0 to ldim-1 do
      for j:=0 to odim-1 do
         for k:=0 to bdim-1 do
             if Criterion[i,j,k]<cmin then begin
                cmin:=Criterion[i,j,k];
                imin:=i;
                jmin:=j;
                kmin:=k;  
             end;
   Lambda:=Grid_LambdaMin+imin*Grid_LambdaStep;
   Omega:=Grid_OmegaMin+jmin*Grid_OmegaStep;
   B:=Grid_BMin+kmin*Grid_BStep;
   Result:=cmin;
   WriteLn('r Minimization results'); 
   WriteLn('r ====================');  
   WriteLn('r Criterion = ',cmin:12:6);
   WriteLn('r Lambda    = ',Lambda:12:4);
   WriteLn('r Omega     = ',Omega:12:8);
   WriteLn('r B         = ',B:12:4); 
   WriteLn('r ====================');
end;

procedure Paint_Slice(Name: String);

var
   TBits: TBitmap;
   j,k,Ampl: Integer;
   SlMin,SlMax,ACoef,BCoef: Double;

begin
   TBits:=TBitmap.Create;
   TBits.PixelFormat:=pf32Bit;
   TBits.Width:=odim;
   TBits.Height:=bdim;
   SlMin:=Slice[0,0];
   SlMax:=SlMin;
   for j:=0 to odim-1 do
      for k:=0 to bdim-1 do
         if Slice[j,k]<SlMin then
            SlMin:=Slice[j,k]
         else if Slice[j,k]>SlMax then
            SlMax:=Slice[j,k];
   ACoef:=255.0/(SlMax-SlMin);
   BCoef:=-ACoef*SlMin;
   with TBits.Canvas do begin
      for j:=0 to odim-1 do
         for k:=0 to bdim-1 do begin
            Ampl:=Round(ACoef*Slice[j,bdim-k-1]+BCoef);
            Ampl:=(Ampl shl 16) + (Ampl shl 8) +Ampl;
            TBits.Canvas.Pixels[j,k]:=TColor(Ampl);
         end
   end;
   TBits.SaveToFile(Name+'.bmp');
   TBits.Destroy;
end;

procedure Compute_Omega_Vector;

var
   j,k: Integer;
   cmi: Double;

begin
   Compute_Omega_B_Slice;
   SetLength(Vector,odim);
   for j:=0 to odim-1 do begin
      cmi:=Slice[j,0];
      for k:=1 to bdim-1 do
         if Slice[j,k]<cmi then
             cmi:=Slice[j,k];
      Vector[j]:=cmi;
   end;
   Vector_Min:=Grid_OmegaMin;
   Vector_Step:=Grid_OmegaStep;
end;

procedure Compute_B_Vector;

var
   j,k: Integer;
   cmi: Double;

begin
   Compute_Omega_B_Slice;
   SetLength(Vector,bdim);
   for k:=0 to bdim-1 do begin
      cmi:=Slice[0,k];
      for j:=1 to odim-1 do
         if Slice[j,k]<cmi then
             cmi:=Slice[j,k];
      Vector[k]:=cmi;
   end;
   Vector_Min:=Grid_BMin;
   Vector_Step:=Grid_BStep;
end;

procedure Compute_Lambda_Vector;

var
   i,j,k: Integer;
   cmi: Double;

begin
   SetLength(Vector,ldim);
   for i:=0 to ldim-1 do begin
      cmi:=Criterion[i,0,0];
      for j:=0 to odim-1 do
         for k:=0 to bdim-1 do
             if Criterion[i,j,k]<cmi then
                cmi:=Criterion[i,j,k];
      Vector[i]:=cmi; 
   end;
   Vector_Min:=Grid_LambdaMin;
   Vector_Step:=Grid_LambdaStep;
end;

procedure Write_Vector(Name: String);

var
   i: Integer;
   arg: Double;
   F: Text;

begin
   Assign(F,Name+'.dat');
   Rewrite(F);
   WriteLn(F,'c Vector');
   WriteLn(F,'h 100000.0');
   WriteLn(F,'l -100000.0');
   for i:=0 to High(Vector) do begin
      arg:=Vector_Min+i*Vector_Step;
      WriteLn(F,arg:12:6,' ',Vector[i]:12:8); 
   end;      
   Close(F);
end;

procedure Write_Profile(Name: String);

var
   i: Integer;
   arg: Double;
   F: Text;

begin
   Assign(F,Name+'.dat');
   Rewrite(F);
   WriteLn(F,'c Profile');
   for i:=0 to 360 do begin
      arg:=i;
      WriteLn(F,arg:12:6,' ',(EvaluatorProfile(arg)):12:6);  
   end;
   Close(F);
end;

procedure Write_Phases(Name: String);

var
   i,j: Integer;
   F: Text;
   Ph,Ampl: Double;

begin
   Assign(F,Name+'.dat');
   Rewrite(F);
   WriteLn(F,'c Phase Ampl Long');
   for i:=rmin to rmax do
      if CR[i] then begin
         for j:=1 to Cnt[i] do begin
            Ph:=Frac((Longitude[i,j]-Shift[i])/360.0);
            if Ph<0.0 then
               Ph:=Ph+1.0;
            Ampl:=Area[i,j]; 
            Ph:=Ph*360.0;
            WriteLn(F,i:4,' ',Ph:10:5,' ',Ampl:10:5,' ',(Longitude[i,j]):10:5); 
         end;
      end; 
   Close(F);
end;

procedure Eliminate_Longitudes;

var
   i,j: Integer;

begin
   for i:=rmin to rmax do begin
      if CR[i] then
         for j:=1 to Cnt[i] do 
             Longitude[i,j]:=0.0;
      Shift[i]:=-Shift[i];
   end;
end;

procedure Write_Distribution(Name: String);

var
   i,j,iadr: Integer;
   Ph,Grand,Cn: Double;
   Distribution: array [0..35] of Double;
   F: Text;

begin
   Assign(F,Name+'.dat');
   Rewrite(F);
   for i:=0 to 35 do
      Distribution[i]:=0.0;
   Grand:=0.0;
   for i:=rmin to rmax do
      if CR[i] then begin
         for j:=1 to Cnt[i] do begin
            Ph:=Frac((Longitude[i,j]-Shift[i])/360.0);
            if Ph<0.0 then
               Ph:=Ph+1.0;
            iadr:=Trunc(36*Ph);
            if iadr>35 then
               iadr:=35
            else if iadr<0 then
               iadr:=0;
            Distribution[iadr]:=Distribution[iadr]+Area[i,j];
            Grand:=Grand+Area[i,j];
         end;     
      end;
   WriteLn(F,'r Total Area = ',Grand:12:4);
   WriteLn(F,'c Dist');
   for i:=0 to 35 do begin
      Cn:=Distribution[i]; 
      WriteLn(F,(10*i):5,' ',Cn:8:4);
   end;
   Close(F);
end;

procedure Compile_Full_Dataset(MinR,MaxR: Integer; IsNorth: Boolean; Name: String);

begin
   Min_Rot:=MinR;
   Max_Rot:=MaxR;
   ReadDataFromSource('all.org',Name+'.dat',IsNorth);
   ReadMeansFromData(Name+'.dat');
   InterpolateMeans;
   WriteMeans(Name+'corr');
   RestoreDefaults;
end;

procedure Compile_Strip(MinR,MaxR: Integer; IsNorth: Boolean; LMin,LMax: Double; Name: String);

begin
   Min_Rot:=MinR;
   Max_Rot:=MaxR;
   Min_Lat:=LMin;
   Max_Lat:=LMax;
   ReadDataFromSource('all.org',Name+'.dat',IsNorth);
   ReadMeansFromData(Name+'.dat');
   InterpolateMeans;
   WriteMeans(Name+'corr');
   RestoreDefaults;
end;

procedure Compute_Method(LambdaMin,LambdaMax,LambdaStep,
                         OmegaMin,OmegaMax,OmegaStep,
                         BMin,BMax,BStep: Double; 
                         Name: String; 
                         Method: TMethod; PhysMode,Repar: Boolean; 
                         var Lambda,Omega,B: Double);

 
begin
   WriteLn('r ================================================================================');
   WriteLn('r Compute Method');
   WriteLn('r');
   WriteLn('r Input  = ',Name);
   ReadMeans(Name+'corr.dat');
   if BMin=BMax then
      Repar:=False; 
   ComputeCorrections;
   if Repar then
      ReparametrizeCorrections;
   ReadData(Name+'.dat',True);
   CheckData(True); 
   case Method of
   mtHits: begin
          WriteLn('r Hit Method');
          Evaluator:=Compute_Hits;
          Norm_Coef:=CountRecords;
      end;
   mtSimple: begin 
          WriteLn('r Simple Method');
          Evaluator:=Compute_Simple;
      end;
   mtDouble: begin
          WriteLn('r Double Method');
          Evaluator:=Compute_Double; 
   end;
   mtShifts: begin
          WriteLn('r Compute From Shifts Method');
          Evaluator:=Compute_From_Shifts; 
   end;
   end;
   Compute_Grid(LambdaMin,LambdaMax,LambdaStep,OmegaMin,OmegaMax,OmegaStep,BMin,BMax,BStep,PhysMode);
   Compute_Grid_Minimum(Lambda,Omega,B);
   if Repar and (OmegaMin<>OmegaMax) then begin
      ReadMeans(Name+'corr.dat');
      ComputeCorrections;
      Compute_Grid(LambdaMin,LambdaMax,LambdaStep,OmegaMin,OmegaMax,OmegaStep,B,B,0.0,PhysMode);  
      Compute_Grid_Minimum(Lambda,Omega,B);   
   end;
   Evaluator:=Compute_Simple;
end;

procedure Write_Shifts(Lambda,Omega,B: Double; Name,ResName: String; PhysMode: Boolean);

begin
   WriteLn('r Write Shifts');
   WriteLn('r Name   = ',Name);
   WriteLn('r Output = ',ResName);
   WriteLn('r Lambda = ',Lambda:12:6);
   WriteLn('r Omega  = ',Omega:12:6);
   WriteLn('r B      = ',B:12:6);
   ReadMeans(Name+'corr.dat');
   ComputeCorrections;
   ComputeShifts(Lambda,Omega,B,PhysMode);
   if ResName<>'' then
      WriteShifts(ResName); 
end;

procedure Write_Corrections(Name,ResName: String);

begin 
   WriteLn('r Write Corrections');
   WriteLn('r Name   = ',Name);
   WriteLn('r Output = ',ResName);
   ReadMeans(Name+'corr.dat');
   ComputeCorrections;
   WriteCorrections(ResName);
end;

procedure Write_Results(Lambda,Omega,B: Double; Name,ResName: String);

var
   NRot,NRec: Integer;
   NormCoef,Cri,NAxi,ACoef,ROmega,Df,Ph: Double;
   
begin
   ReadMeans(Name+'corr.dat');
   ComputeCorrections;
   ACoef:=SlopeOfCorrections;
   ReadData(Name+'.dat',True);
   CheckData(True); 
   NormCoef:=ComputeNorm;   
   NRot:=CountRotations;
   NRec:=CountRecords;
   ComputeShifts(Lambda,Omega,B,False);
   Cri:=Compute_Simple/NormCoef;
   NAxi:=Compute_NAxi;
   if ResName<>'' then begin
      Write_Phases(ResName);
      Write_Distribution(ResName+'D');
   end;
   WriteLn;
   Write('Results for data set ',Name,'.dat');
   if ResName<>'' then
       WriteLn(', phases in ',ResName,'.dat, distribution in ',ResName+'D.dat.')
   else
       WriteLn;
   WriteLn;
   WriteLn('Omega  = ',Omega:12:6);
   WriteLn('B      = ',B:12:4);
   WriteLn('Lambda = ',Lambda:12:2); 
   WriteLn;
   ROmega:=Omega-B*ACoef;
   Df:=360.0/(360.0/25.38-ROmega);
   Df:=Df/365.25;
   Ph:=2*360.0/25.38-Omega;
   WriteLn('Mean rotation rate        = ',ROmega:12:6);
   WriteLn('Number of extra rotations = ',((Abs(Shift[rmax]-Shift[rmin]))/360.0):12:4);
   WriteLn('Cycle Length              = ',Df:12:6);
   WriteLn('Physical Omega            = ',Ph:12:6); 
   WriteLn('Physical B                = ',(-B):12:6); 
   WriteLn('Slope from waves          = ',(B*ACoef):12:6);
   WriteLn;
   WriteLn('Criterion       = ',Cri:12:6);
   WriteLn('Non-Axisymmetry = ',NAxi:12:6);
   WriteLn;

   WriteLn('Number of used rotations = ',NRot:5);
   WriteLn('Number of used records   = ',NRec:6); 
   WriteLn;
   WriteLn('******************************************************************************************************************************');
end;

procedure RestoreDefaults;

begin
   Is_Corrected:=True;
   Include_Singles:=True;
   Min_Rot:=1;
   Max_Rot:=Max_Rotations;
   Min_Year:=0;
   Max_Year:=2500; 
   JD_Min:=0.0;
   JD_Max:=3000000.0;
   Min_Lat:=0.0;
   Max_Lat:=90.0;
   Min_Area:=0;
   Max_Area:=100000; 
end;

function Compute_Hits: Double;

var
   i,j,Sum,Hit_R: Integer;
   a,b,Ph,D: Double;
   Hits: array of array of Boolean;

begin
   Hit_R:=Trunc((rmax-rmin)/Rot_Count)+1;
   a:=(Hit_R-1.0)/(rmax-rmin);
   b:=-a*rmin;
   SetLength(Hits,Hit_R,Hit_D);   
   for i:=0 to Hit_R-1 do
      for j:=0 to Hit_D-1 do
         Hits[i,j]:=False;    
   D:=Hit_D-1;
   for i:=rmin to rmax do
      if CR[i] then begin
         for j:=1 to Cnt[i] do begin
            Ph:=Frac((Longitude[i,j]-Shift[i])/360.0);
            if Ph<0.0 then
               Ph:=Ph+1.0;
            Hits[Round(a*i+b),Round(D*Ph)]:=True;
         end;
      end;
    Sum:=0;
    for i:=0 to Hit_R-1 do
       for j:=0 to Hit_D-1 do
           if Hits[i,j] then
               Inc(Sum);
    Result:=Sum/Norm_Coef;      
end;

function Compute_Simple: Double;

var
   i,j: Integer;
   Ph: Double;
   
begin
   Result:=0.0;
   for i:=rmin to rmax do
      if CR[i] then begin
         for j:=1 to Cnt[i] do begin
            Ph:=Frac((Longitude[i,j]-Shift[i])/360.0);
            if Ph<0.0 then
               Ph:=Ph+1.0;
            if Ph<0.5 then begin
               if Ph<0.25 then
                  Ph:=-Ph+0.25
               else
                  Ph:=Ph-0.25;  
            end else begin
               if Ph<0.75 then
                  Ph:=-Ph+0.75
               else
                  Ph:=Ph-0.75;   
            end;
            Result:=Result+Area[i,j]*Sqr(Ph);
         end;
     end;
end;

function Compute_Simple_Profile(Deg: Double): Double;

var
   Ph: Double;

begin
   Ph:=Frac(Deg/360);
   if Ph<0.0 then
      Ph:=Ph+1.0; 
   if Ph<0.5 then begin
      if Ph<0.25 then
         Ph:=-Ph+0.25
      else
         Ph:=Ph-0.25;  
   end else begin
      if Ph<0.75 then
         Ph:=-Ph+0.75
      else
         Ph:=Ph-0.75;   
   end;
   Result:=Ph;
end;

function Compute_From_Shifts: Double;

var
   i: Integer;
   Ph: Double;
   
begin
   Result:=0.0;
   for i:=rmin to rmax do
      if CR[i] then begin 
         Ph:=Frac((-Shift[i])/360);
         if Ph<0.0 then
            Ph:=Ph+1.0; 
         if Ph<0.5 then begin
            if Ph<0.25 then
               Ph:=-Ph+0.25
            else
               Ph:=Ph-0.25;  
         end else begin
            if Ph<0.75 then
               Ph:=-Ph+0.75
            else
               Ph:=Ph-0.75;   
         end;
         Result:=Result+Sqr(Ph);
      end;
end;

function Compute_NAxi: Double;

var
   i,j: Integer;
   Ph,N1,N2: Double;
   
begin
   N1:=0.0;
   N2:=0.0;
   for i:=rmin to rmax do
      if CR[i] then begin
         for j:=1 to Cnt[i] do begin
            Ph:=Frac((Longitude[i,j]-Shift[i])/360);
            if Ph<0.0 then
               Ph:=Ph+1.0; 
            if (Abs(Ph-0.25)<0.125) or (Abs(Ph-0.75)<0.125) then
               N1:=N1+Area[i,j]
            else
               N2:=N2+Area[i,j];
         end;
   end;
   Result:=(N1-N2)/(N1+N2);
end;

function Compute_Double: Double;

var
   i,j: Integer;
   Ph: Double;
   
begin
   Result:=0.0;
   for i:=rmin to rmax do
      if CR[i] then begin
         for j:=1 to Cnt[i] do begin
            Ph:=Frac((Longitude[i,j]-Shift[i])/360);
            if Ph<0.0 then
               Ph:=Ph+1.0;
            Ph:=Frac(Ph*4.0); 
            Ph:=Abs(0.5-Ph);
            Result:=Result+Area[i,j]*Sqr(Ph);
         end;
   end;
end;

initialization
   RestoreDefaults;
   ldim:=0;
   odim:=0;
   bdim:=0;
   Rot_Count:=10;
   Hit_D:=60;
   Norm_Coef:=1.0;
   Evaluator:=Compute_Simple;
   EvaluatorProfile:=Compute_Simple_Profile;
end.