unit CloudRemUnit;

{$mode objfpc}{$H+}

interface

uses
  Classes, SysUtils, LResources, Forms, Controls, Graphics, Dialogs, StdCtrls,
  Buttons, ExtCtrls, ComCtrls
  , strutils
  , DateUtils
  , LazFileUtils //required to extractfilenameonly for converting to csv files
  , moon //required for Moon calculations
  , Math
  , header_utils;

type

  { TCloudRemMilkyWay }

  TCloudRemMilkyWay = class(TForm)
    Label2: TLabel;
    HalfRangeEdit: TLabeledEdit;
    LocationEdit: TLabeledEdit;
    LatEdit: TLabeledEdit;
    LongEdit: TLabeledEdit;
    Memo1: TMemo;
    SourceFileDialog: TOpenDialog;
    ProcessStatusMemo: TMemo;
    ProgressBar1: TProgressBar;
    SourceFileButton: TBitBtn;
    SourceFileEdit: TEdit;
    StartButton: TButton;
    StatusBar1: TStatusBar;
    procedure FormShow(Sender: TObject);
    procedure SourceFileButtonClick(Sender: TObject);
    procedure StartButtonClick(Sender: TObject);
  private

  public
    procedure FileTypeCheck();
  end;

var
  CloudRemMilkyWay: TCloudRemMilkyWay;
  SourceFileName: String;
  Filetype: String = 'unknown';
  SQM_Lat, SQM_Long: Extended;
  LocationName: String;
  half_range: integer;
  OutputHeader: String;

implementation
uses
  appsettings, Unit1
;
{$R *.lfm}

Const
Section='CloudRemovalMilkyWayTool';

procedure Memo(s:String);
begin
  CloudRemMilkyWay.ProcessStatusMemo.Append(s);
  CloudRemMilkyWay.ProcessStatusMemo.SelStart:=CloudRemMilkyWay.ProcessStatusMemo.GetTextLen;
  Application.ProcessMessages;
end;

procedure TCloudRemMilkyWay.FileTypeCheck();
type
   TFloatArray = Array [1..100] of Float;
var
     Str: String;
     pieces: TStringList;
     fdata: TextFile;
     PositionValid: Boolean = True; //Assume valid location
     RequiredCSVFields: Integer = 21;
     LineCount: Integer = 0;
     FileCheckValid: Boolean = True; //Assume valid file data
     MeanArray : TFloatArray;
     OldUTCRecord, NewUTCRecord :TDateTime;
begin
  pieces := TStringList.Create;

  if AnsiContainsStr(SourceFileName, '.csv') then Filetype:='csv'
  else if AnsiContainsStr(SourceFileName, '.dat') then Filetype:='dat'
  else Filetype:='unknown';

  case Filetype of
         'dat': begin
           LatEdit.Enabled:=False;
           LongEdit.Enabled:=False;

           AssignFile(fdata, SourceFileName);
           reset(fdata);

           //Read first line
           Readln(fdata, Str);

           while AnsiStartsStr('#',Str) do begin
             // Get location data from header.

             if AnsiStartsStr('# Location name:',Str) then begin
                 //Remove comment from beginning of line.
                 LocationName := AnsiRightStr(Str,length(Str) - RPos(':',Str));

                 //Remove any text after comma
                 LocationName := AnsiLeftStr(LocationName,length(LocationName) - RPos(',',LocationName));

                 //Remove surrounding spaces and Replace spaces with underscore
                 LocationName:=StringReplace(Trim(LocationName),' ','_',[rfReplaceAll]);

                 if (Length(LocationName)=0) then
                     LocationName:= 'Not-Specified';

                 LocationEdit.Text:=LocationName;
                 LocationEdit.Enabled:=False;

             end;


             if AnsiStartsStr('# Position',Str) then begin
                 //Prepare for parsing.
                 pieces.Delimiter := ',';
                 pieces.StrictDelimiter := True; //Do not parse spaces.
                 //Remove comment from beginning of line.
                 pieces.DelimitedText := AnsiRightStr(Str,length(Str) - RPos(':',Str));

                 //Check for location entered
                 if (pieces.Count<2) then begin
                   Memo('Error: No position data found in header, an example is:');
                   Memo('   # Position (lat, lon, elev(m)): 43.24611, -118.8942, 1256');
                   LatEdit.Text:='nul';
                   LongEdit.Text:='nul';
                   PositionValid:=False;
                   SQM_Lat:=0;
                   SQM_Long:=0;
                   StartButton.Enabled:=False;
                 end
                 else
                   begin
                     StartButton.Enabled:=True;
                     SQM_Lat:=StrToFloat(pieces.Strings[0]);
                     if ((SQM_Lat<-90.0) or (SQM_Lat>90.0)) then begin
                        Memo('Error: Latitude ' + FloatToStr(MyLatitude) +' out of range (-90 to 90).');
                        PositionValid:=False;
                     end;
                     SQM_Long:=StrToFloat(pieces.Strings[1]);
                     if ((SQM_Long<-180.0) or (SQM_Long>180.0)) then begin
                       Memo('Error: Longitude '+FloatToStr(MyLongitude)+' out of range (-180 to 180)');
                       PositionValid:=False;
                     end;
                     LatEdit.Text:=Format('%.7f',[SQM_Lat]);
                     LongEdit.Text:=Format('%.7f',[SQM_Long]);
                   end;

                   StartButton.Enabled:=PositionValid;

               end;

             //Read next line
             Readln(fdata, Str);

           end;


           //Read next 100 lines to determine sample rate
           Memo('Reading first 100 lines to determine sampling rate.');
           Application.ProcessMessages;
           try
             pieces.Delimiter := ';';

             //Throw away first line since it may have started recording on power up instead of at timed interval.
             Readln(fdata, Str);
             pieces.DelimitedText :=Str;
             //Parse out time
             NewUTCRecord:=ScanDateTime('yyyy-mm-dd"T"hh:nn:ss.zzz',pieces.Strings[0]);

             for LineCount:=1 to 100 do begin
               Readln(fdata, Str);
               pieces.DelimitedText :=Str;
               //Move previous new-time to old-time
               OldUTCRecord:=NewUTCRecord;

               //Parse out time
               NewUTCRecord:=ScanDateTime('yyyy-mm-dd"T"hh:nn:ss.zzz',pieces.Strings[0]);
               //Memo(Format('Seconds=%d',[SecondOf(NewUTCRecord)]));
               //Determine time difference, and put into the array for mean calculation.
               //Memo(Format('Seconds difference=%d',[SecondsBetween(NewUTCRecord, OldUTCRecord)]));
               MeanArray[LineCount]:=SecondsBetween(NewUTCRecord, OldUTCRecord);
             end;
             Memo(Format('Mean value = %f second(s) frequency.',[mean(MeanArray)]));
           except
             Memo('Error: Failed to read first 100 lines to determine sampling rate.');
             FileCheckValid:=False;
           end;

           if FileCheckValid then begin
             Memo('Finished reading first 100 lines to determine sampling rate.');
           end;

           StartButton.Enabled:=FileCheckValid;

           CloseFile(fdata);

         end;


         'csv': begin
           AssignFile(fdata, SourceFileName);
           reset(fdata);

           //Read first line
           Readln(fdata, Str);

           //Check for required variables
           //Prepare for parsing.
           pieces.Delimiter := ',';
           pieces.StrictDelimiter := True; //Do not parse spaces.
           //Remove comment from beginning of line.
           pieces.DelimitedText := AnsiRightStr(Str,length(Str) - RPos(':',Str));

           //Check for proper field count
           if pieces.Count=RequiredCSVFields then begin //Proper number of fields in.csv file
             LocationEdit.Enabled:=True;

             //retrieve location name from udm config file
             LocationName:=vConfigurations.ReadString(Section,'LocatioName','0');
             LocationEdit.Text:=LocationName;

             SQM_Lat:= StrToFloat(vConfigurations.ReadString(Section,'Latitude','0'));
             LatEdit.Text:=Format('%.7f',[SQM_Lat]);

             SQM_Long:=StrToFloat(vConfigurations.ReadString(Section,'Longitude','0'));
             LongEdit.Text:=Format('%.7f',[SQM_Long]);

             half_range:=StrToInt(vConfigurations.ReadString(Section,'half_range','9'));
             HalfRangeEdit.Text:=Format('%d',[half_range]);
           end
           else begin  //Improper number of fields in.csv file
             Memo(Format('Error: Improper number of fields %d in.csv file, %d required.',[pieces.Count, RequiredCSVFields]));

             Memo('Got:');
             Memo(Str);
             Memo('');
             Memo('Should be:');
             Memo('Location,UTC_YYYY,UTC_MM,UTC_DD,UTC_HH,UTC_mm,UTC_ss,Local_YYYY,Local_MM,Local_DD,Local_HH,Local_mm,Local_ss,Celsius,Volts,Msas,Status,MoonPhaseDeg,MoonElevDeg,MoonIllum%,SunElevDeg');

             LatEdit.Enabled:=False;
             LongEdit.Enabled:=False;

             LatEdit.Text:='nul';
             LongEdit.Text:='nul';
             PositionValid:=False;
             SQM_Lat:=0;
             SQM_Long:=0;
             StartButton.Enabled:=False;
           end;
           StartButton.Enabled:=PositionValid;
           CloseFile(fdata);
         end;

         else begin
           LocationEdit.Enabled:=True;
           LatEdit.Enabled:=True;
           LongEdit.Enabled:=True;
         end;
  end;

end;

  function yisleap(year:integer):Integer;
  begin
    if  (year mod 4=0)
    and ((year mod 100>0) or (year mod 400 = 0)) then result:=1
   else result:=0;
  end;

  function get_yday(mon: integer;  day: integer;  year: integer): integer;
  const
  days: array [0..Pred(2),0..Pred(13)] of integer = ((0,0,31,59,90,120,151,181,212,243,273,304,334),(0,0,31,60,91,121,152,182,213,244,274,305,335));
  var
  days_this_year: integer;
  n: integer;
  days_since_Jan1_2018: integer;
  sum: integer;

  leap: integer;

  begin
    sum:=0;
    leap:=yisleap(year);
    days_this_year:= days[leap][mon]+day; (* count the number of days in complete previous years *)
    (* note, this algorithm does not work for any dates prior to 1-1-2018 *)
    for n:= year-1 downto 2018 do
    begin
      leap:=yisleap(n);
      if leap=0 then begin
        sum:= sum+365;
      end
      else begin
        sum:= sum+366;
      end;
    end;
    days_since_Jan1_2018:= sum+days_this_year;
    //Memo(Format(' Year=%d Days in previous year(s)=%d',[year,sum]));
    //Memo(Format(' Month=%d Date=%d Days this year=%d',[mon, day, days_this_year]));
    begin
      result:= days_since_Jan1_2018;
      exit;
    end;
  end;

  function get_UT(UTC_Hour: integer;  UTC_Min: integer;  UTC_Sec: integer): double;
  var
     UT: double;
  begin
    UT:= UTC_Hour+ double(UTC_Min) / 60.0 + double(UTC_Sec) / 3600.0;
    //Memo(Format('UTC hour=%d, UTC_Min=%d, UTC_Sec=%d',[UTC_Hour, UTC_Min, UTC_Sec]));
    //Memo(Format(' in_get_UT  UT= %7.6f',[UT]));
    begin
      result:= UT;
      exit;
    end;
  end;

  function get_J2000(year: integer;  UTC_Month: integer;  UTC_Day: integer;  UTC_Hour: integer;  UTC_Min: integer;  UTC_Sec: integer): double;
  var
  dwhole: double;
  dfrac: double;
  J2000_days: double;
  begin
      dwhole := 367 * year - (7*(year + (UTC_Month + 9) div 12) div 4) + (275 * UTC_Month div 9) + UTC_Day - 730531.5;
      dfrac := (double(UTC_Hour) + (double(UTC_Min))/60. + (double(UTC_Sec))/3600.0)/24.0;
      J2000_days := dwhole + dfrac;
      //Memo(Format(' in get_J2000  J2000_days=%10.6f',[J2000_days]));
      Result := J2000_days;
  end;



  function get_right_ascension(year: integer;  UTC_Month: integer;  UTC_Day: integer;  UTC_Hour: integer;  UTC_Min: integer;  UTC_Sec: integer;  SQM_Long: double): double;
  var
  right_ascension: double;
  J2000_days: double;
  UT: double;
  //aa: double;
  //bb: double;
  //cc: double;
  //dd: double;
  //ee: double;
  multiples: integer; (* need to calculate Universal Time (UT) and J2000_day *)
  (* set up test data for which we know the answer
      year = 2008;
      UTC_Month = 4;
      UTC_Day = 4;
      UTC_Hour = 15;
      UTC_Min = 30;
      UTC_Sec = 0.0;
      SQM_Long = -1.9166667;
  *)
  (* get the J2000 day value *)
  begin

    J2000_days:= get_J2000(year,UTC_Month,UTC_Day,UTC_Hour,UTC_Min,UTC_Sec);
    //Memo(Format(' in get_right_ascension -- SQM_Long= %6.3f',[SQM_Long]));
    UT:= get_UT(UTC_Hour,UTC_Min,UTC_Sec);
    //Memo(Format(' J2000_days= %6.2f',[J2000_days]));
    (*  get the Universal time as a fraction of a day *)
    //Memo(Format(' UT= %6.4f',[UT]));

    right_ascension:= 100.46+0.985647*J2000_days+SQM_Long+15.*UT;
    //Memo(Format(' right_ascension before 0-360 check(degrees) = %6.2f',[right_ascension]));
    multiples:= Round(right_ascension / 360.0); (* make sure that the value is within the range of 0 to 360 degrees *)
    (* how many multiples of 360 do we have?  Subtract or add out that number *)
    (* note that multiples is an integer and that the remainder is truncated in the next calculation *)
    if multiples>0 then begin
      right_ascension:= right_ascension- double(multiples)*360.0;
    end;
    if multiples<0 then begin
      right_ascension:= right_ascension- double(multiples)*360.0;
    end;
    if right_ascension<0 then begin
      right_ascension:= right_ascension+360.;
    end;
    //Memo(Format(' right_ascension after 0-360 check(degrees) = %6.3f',[right_ascension]));
    right_ascension:= right_ascension / 15;
    //Memo(Format(' right_ascension (hours) = %6.4f',[right_ascension]));
    (* convert right_ascension from degrees to hours *)
    begin
      result:= right_ascension;
      exit;
    end;
  end;


  //function addSQMattr(argc : integer;var argv : byte):integer;
  procedure addSQMattr();
  begin
  end;

  { TCloudRemMilkyWay }

  procedure TCloudRemMilkyWay.StartButtonClick(Sender: TObject);
    var
    k,m: integer;

    minutes_since_3pm: array [0..Pred(1500)] of integer;
    dUYear: array [0..Pred(1500)] of integer;
    dUMonth: array [0..Pred(1500)] of integer;
    dUDay: array [0..Pred(1500)] of integer;
    dUHour: array [0..Pred(1500)] of integer;
    dUMinute: array [0..Pred(1500)] of integer;
    dUSeconds: array [0..Pred(1500)] of double;
    dYear: array [0..Pred(1500)] of integer;
    dMonth: array [0..Pred(1500)] of integer;
    dDay: array [0..Pred(1500)] of integer;
    dHour: array [0..Pred(1500)] of integer;
    dMinute: array [0..Pred(1500)] of integer;
    dSeconds: array [0..Pred(1500)] of double;
    dMsas: array [0..Pred(1500)] of double;
    //dMsas_Corr: array [0..Pred(1500)] of double;
    dVolts: array [0..Pred(1500)] of double;
    dCelsius: array [0..Pred(1500)] of double;
    dMoonPhase: array [0..Pred(1500)] of double;
    dMoonElev: array [0..Pred(1500)] of double;
    dMoonIllum: array [0..Pred(1500)] of double;
    dSunElev: array [0..Pred(1500)] of double;
    msas_Avg: array [0..Pred(1500)] of double;
    msas_Sum: double;
    msas_Count: double;
    dStatus: array [0..Pred(1500)] of integer;
    //NameIn: array [0..Pred(120)] of char;
    NameOut: array [0..Pred(250)] of char;
    //SQM_Location: array [0..Pred(30)] of char;
    //blank: array [0..Pred(200)] of char;
    //nfile: integer;
    //Slength: integer;
    //ret: integer;
    Start: integer;
    Last: integer;
    len2: integer;
    right_ascension: double;
    SQM_Dec, SQM_RA: Extended;
    J2000_days: double;
    timediff: integer;
    num_minutesA: integer;
    num_minutesB: integer;
    //remainder: double;
    Galactic_Lat: double;
    Galactic_Long: double;
    XX, YY: double;
    //Galactic_Elevation: double;
    pi: double;

    (* NGP is North Galactic Pole, NCP is North Celestial Pole *)
    RightAscension_NGP, Dec_NGP, Galactic_Long_NCP: double;
    sum_x, sum_y, sum_xy, sum_x2, sum_y2: Extended;
    N, mean_x, mean_y, mean_xy, mean_x2: Extended;
    slope, yintercept, rcorr, rsqrd: Extended;
    kk, timediff_max: integer;
    Observed: array [0..Pred(1500)] of Extended;
    Expected: array [0..Pred(1500)] of Extended;
    nodata1: Extended;
    nodata2: Extended;
    DOF: Extended;
    RSE: array [0..Pred(1500)] of Extended;
    SS: array [0..Pred(1500)] of Extended;
    RSE_mult: Extended; (* set up a multiplier for the RSE values and no-data output values *)
    (* we output the RSE values multiplied by this constant to give more manageable values *)
    fdata, fdataout: TextFile;
    outstr: String;

    (* added to handle the daylight savings time fix to "minutes since 3pm" *)
    dPosNeg, dHour_Delta, dShift_Hour: Integer;

    days: integer;
     (* We actually want the number of nights since Jan 1, 2018 - that is we want to count the evening and night as part of the
     * same "day" - actually the same "night"; So if the minutes since 15:00 hours is greater than 540 (i.e. after midnight) we
     * subtract one day from the "days" value so those times are considered part of the previous day (i.e."night") *)

     Str: String;
     pieces: TStringList;

     //.dat fields, -1 = not defined  yet.
     UTCTimeField: Integer =-1; //Field that contains the UTC time variable.
     LocalTimeField: Integer =-1; //Field that contains the Local time variable.
     TemperatureField: Integer = -1; //Field that contains the Temperature variable.
     VoltageField: Integer = -1; //Field that contains the Voltage variable.
     MSASField: Integer = -1; //Field that contains the MSAS variable.
     RecordTypeField: Integer = -1; //Field that contains the Record Type (Initial/subsequent) variable.

     SourceExtension: String; //.dat or .csv

     UTCRecord :TDateTime;
     LocalRecord :TDateTime;

     MoonElevation: extended = 0.0;
     MoonAzimuth: extended = 0.0;
     SunElevation: extended = 0.0;
     SunAzimuth: extended = 0.0;

     i: Integer = 0; //General purpose counter
     ErrorInputLineCounter:Int64;
     ExceptionFlag: Boolean = False;



    label ReadAnother, LastDay, Termination;

    begin

      ProcessStatusMemo.Clear;

      pieces := TStringList.Create;
      pieces.Delimiter := ',';
      //pieces.StrictDelimiter := False; //Parse spaces also
      pieces.StrictDelimiter := True; //Do not parse spaces.

      k:=0;
      m:=0;
      ErrorInputLineCounter:=0;

      RSE_mult:= 1000.;
      nodata1:= 999.*RSE_mult;
      nodata2:= 888.*RSE_mult;
      timediff_max:= 16;
      Memo(Format('We are running Program %s',[ParamStr(0)]));
      (* specify the maxiumum number of minutes allowed between SQM samples, prior to marking a data gap *)
      (* Run this program by specifying the program name, followed by three parameters:
       *                      1) A file of SQM data which has already been processed as a csv with sun and moon data
       *                      2) The latitude of the SQM location in fractions and
       *                      3) The longitude of the SQM location in fractions and so the command line should look like this:
       *                             ./addSQMattributes inputfilename.csv 43.7916667 -120.23422 *)



      len2:= length(SourceFileName);
      //Memo(Format(' The input filename on reading is: %s which has %d characters.',[SourceFileName,len2-1]));
      Memo(Format(' The input filename is: %s',[SourceFileName]));


      //Save LocationName
      vConfigurations.WriteString(Section,'LocatioName',LocationName);

      //Get Latitude
      SQM_Lat:=StrToFloat(LatEdit.Text);
      Memo(Format(' The latitude of the SQM is: %.7f',[SQM_Lat]));
      //Save Latitude
      vConfigurations.WriteString(Section,'Latitude',Format('%.7f',[SQM_Lat]));

      //Get Longitude
      SQM_Long:=StrToFloat(LongEdit.Text);
      Memo(Format(' The longitude of the SQM is: %.7f',[SQM_Long]));
      //Save Longitude
      vConfigurations.WriteString(Section,'Longitude',Format('%.7f',[SQM_Long]));

      half_range:=StrToIntDef(HalfRangeEdit.Text,0);
      Memo(Format(' The Half Range is: %d',[half_range]));
      vConfigurations.WriteString(Section,'half_range',Format('%d',[half_range]));

      (* Open the input file *)
      AssignFile(fdata, SourceFileName);
      reset(fdata);

      SourceExtension:=ExtractFileExt(SourceFileName);

      (* Open an output file to hold the output data *)
      (* tack on "SQM_attr" before the .csv *)
      NameOut:=concat( ExtractFileNameWithoutExt(SourceFileName),'_sun-moon-mw-clouds.csv');
      //Memo(concat( ExtractFileNameWithoutExt(SourceFileName),'_sun-moon-mw-clouds.csv'));
      Memo(Format('The Output Data Filename is: %s',[NameOut]));
      AssignFile(fdataout, NameOut);
      Rewrite(fdataout);//Create the file


      case SourceExtension of
             '.dat': begin
               (*  Get the field locations from the .dat file
                   The bare minimum is:
                      UTC Date & Time, Local Date & Time, Temperature, Voltage, MSAS, Record type
                   The desired output is:
                      Location,UTC_YYYY,UTC_MM,UTC_DD,UTC_HH,UTC_mm,UTC_ss,Local_YYYY,Local_MM,Local_DD,Local_HH,Local_mm,Local_ss,Celsius,Volts,Msas,Status,MoonPhaseDeg,MoonElevDeg,MoonIllum%,SunElevDeg  *)

               Readln(fdata, Str);
               while AnsiStartsStr('#',Str) do begin
                 (* Read the data file *)

                 if AnsiStartsStr('# UTC Date & Time',Str) then begin
                   pieces.Delimiter := ',';
                   pieces.DelimitedText :=Str;
                   //Memo(Format(' Str=%s, i=%d, pieces.Count=%d',[Str, i, pieces.Count]));
                    for i:=0 to pieces.Count-1 do begin
                      case Trim(pieces.Strings[i]) of
                        '# UTC Date & Time': UTCTimeField:=i;
                        'Local Date & Time': LocalTimeField:=i;
                        'Temperature': TemperatureField:=i;
                        'Voltage': VoltageField:=i;
                        'MSAS': MSASField:=i;
                        'Record type': RecordTypeField:=i;
                      end;
                   end;
                    //Memo(Format(' %d, %d, %d, %d, %d, %d ',[UTCTimeField, LocalTimeField, TemperatureField, VoltageField, MSASField, RecordTypeField]));

                    //Read two more useless lines then move to next stage
                    Readln(fdata, Str);
                    Readln(fdata, Str);

                    OutputHeader:='Location,Lat,Long,UTC_Date,UTC_Time,Local_Date,Local_Time,Celsius,';
                    OutputHeader:=OutputHeader+'Volts,';
                    OutputHeader:=OutputHeader+'Msas,';
                    OutputHeader:=OutputHeader+'Status,';
                    OutputHeader:=OutputHeader+'MoonPhase,MoonElev,MoonIllum,SunElev,MinSince3pmStdTime,Msas_Avg,NightsSince_1118,RightAscensionHr,Galactic_Lat,Galactic_Long,J2000days,ResidStdErr';

                    break;
                   end;
                 Readln(fdata, Str);
               end;

             end;

             '.csv': begin

               (* Read the first header record and throw it away *)
               Readln(fdata, Str);

               // Set header string
               OutputHeader:='Location,Lat,Long,UTC_Date,UTC_Time,Local_Date,Local_Time,Celsius,Volts,Msas,Status,MoonPhase,MoonElev,MoonIllum,SunElev,MinSince3pmStdTime,Msas_Avg,NightsSince_1118,RightAscensionHr,Galactic_Lat,Galactic_Long,J2000days,ResidStdErr';
             end


             else begin
               Memo(Format('Cannot deal with extension: %s',[SourceExtension]));
               goto Termination;
             end;
      end;

      N:= Extended((2*half_range)+1);
      Memo('');
      (* half_range is the number of samples (usually 5 minutes apart) to include before and after the
       * current point at which the Residual Standard Error of regression calculation is performed; the full width of the interval in terms of number of points
       * to consider in the calculation is (2*half_range + 1) so that a
      *  half_range of 6 incorporates a total range of 60 minutes;
      *  note that we don't increment by one in the minutes calculation
      *  for the middle point in the time range because it is already taken into account
      *  half_range of 9 evaluates to a 90 minute range
      *)
      Memo(Format(' The half_range parameter is set to: %d',[half_range]));
      Memo(Format(' This means that the Residual Error calculation operates over %d samples',[Round(N)]));
      Memo(Format(' In other words, if the sample spacing is 1 minutes, then the range is %d minutes.',[Round(half_range*2*1)]));
      Memo(Format(' Or              if the sample spacing is 5 minutes, then the range is %d minutes.',[Round(half_range*2*5)]));
      Memo(Format(' Or              if the sample spacing is 15 minutes, then the range is %d minutes.',[Round(half_range*2*15)]));
      Memo('');
      Memo(Format(' Residual Standard Error values that we output are multiplied by %d to achieve larger values.',[Round(RSE_mult)]));
      Memo('');
      Memo(Format(' We allow gaps of %d minutes between SQM samples prior to marking a data gap.',[timediff_max]));
      Memo('');
      (* Write a header record to the output file *)
      writeln(fdataout,OutputHeader);
      pi:= 3.14159265359;
      RightAscension_NGP:= (192.85948*(pi) / 180.0);
      //Memo(Format('RightAscension_NGP=%7.6f',[RightAscension_NGP]));
      (* set up some constant values used later to calculate the Galactic Coordinates of the normal at the SQM location *)
      (* These are from the Wikipedia Celestial coordinate system  *)
      (* convert these constants from degrees to radians *)
      Dec_NGP:= (27.12825*(pi) / 180.0);
      //Memo(Format('Dec_NGP= %7.6f',[Dec_NGP]));
      Galactic_Long_NCP:= 122.93192*(pi) / 180.0;
      //Memo(Format('Galactic_Long_NCP= %7.6f',[Galactic_Long_NCP]));
      Memo('Processing file, please wait ...');


      (* Note that the string read format statement, reads up to the first carriage return in the input file, then reads the carriage return itself *)
      m:= -1;
      Start:= 0; (* initiate the record counter *)
      (* initiate the flag on 15 hundred hour *)

    ReadAnother:
      m:= m+1;   (* increment the counter *)

      //Memo(Format('ReadAnother: m=%d', [m]));

      if m>1499 then begin
        Memo('We have more than 1500 samples for this day.');
        Memo('If this is good data, sorry but this option does not handle a data cadence smaller than 1 minute.');
        Memo('Alternately, does the input file have bad data?');
        Memo('  Was the SQM battery dying and taking a sample too frequently or off schedule?');
        Memo('  The input data are therefore suspect.');
        goto Termination;
      end;

      Readln(fdata, Str);
      //Memo(Format('Str= %s',[Str]));

      if (Length(Str)>0) then begin


        case SourceExtension of
               '.dat': begin
                 pieces.Delimiter := ';';
                 pieces.DelimitedText :=Str;

                 try
                   UTCRecord:=ScanDateTime('yyyy-mm-dd"T"hh:nn:ss.zzz',pieces.Strings[0]);
                   LocalRecord:=ScanDateTime('yyyy-mm-dd"T"hh:nn:ss.zzz',pieces.Strings[1]);
                   dUYear[m]:=YearOf(UTCRecord);
                   dUMonth[m]:=MonthOf(UTCRecord);
                   dUDay[m]:=DayOf(UTCRecord);
                   dUHour[m]:=HourOf(UTCRecord);
                   dUMinute[m]:=MinuteOf(UTCRecord);
                   dUSeconds[m]:=SecondOf(UTCRecord);
                   dYear[m]:=YearOf(LocalRecord);
                   dMonth[m]:=MonthOf(LocalRecord);
                   dDay[m]:=DayOf(LocalRecord);
                   dHour[m]:=HourOf(LocalRecord);
                   dMinute[m]:=MinuteOf(LocalRecord);
                   dSeconds[m]:=SecondOf(LocalRecord);
                 except
                   Memo(Format(' Problem with UTC or Local time: %s',[Str]));
                   ExceptionFlag:=True;
                 end;
                 if ExceptionFlag then goto Termination;

                 try
                   dCelsius[m]:=StrToFloatDef(pieces[TemperatureField],0);
                   if VoltageField>0 then dVolts[m]:=StrToFloat(pieces[VoltageField]);
                 except
                   Memo(Format(' Problem with data in record: %s',[Str]));
                   ExceptionFlag:=True;
                 end;
                 if ExceptionFlag then goto Termination;

                 if (pieces[MSASField]<>'') then
                   dMsas[m]:=StrToFloatDef(pieces[MSASField],0)
                 else begin
                   Inc(ErrorInputLineCounter,1);
                   Memo(Format('We found a bad input record at Local %04d-%.2d-%.2dT%.2d:%.2d:%.2d and are skipping it.',
                     [dYear[m],dMonth[m],dDay[m],dHour[m],dMinute[m],Round(dSeconds[m])]));
                   m:=m-1;
                   goto ReadAnother;
                 end;
                 if RecordTypeField>0 then dStatus[m]:=StrToInt(pieces[RecordTypeField]);
                 dMoonPhase[m]:=moon_phase_angle(StrToDateTime(DateTimeToStr(UTCRecord)));


                 //Calculate Moon position
                 //Change sign for Moon calculations
                 Moon_Position_Horizontal(
                     StrToDateTime(DateTimeToStr(UTCRecord)),
                     -1.0*SQM_Long,
                     SQM_Lat,
                     MoonElevation,
                     MoonAzimuth);

                 dMoonElev[m]:=MoonElevation;
                 dMoonIllum[m]:=current_phase(StrToDateTime(DateTimeToStr(UTCRecord)))*100.0;

                 Sun_Position_Horizontal(
                    StrToDateTime(DateTimeToStr(UTCRecord)),
                    -1.0*SQM_Long,
                    SQM_Lat,SunElevation,
                    SunAzimuth);

                 dSunElev[m]:=SunElevation;


               end;
               '.csv':begin
                 pieces.Delimiter := ',';
                 pieces.DelimitedText :=Str;

                 //SQM_Location:=pieces[0];
                 LocationName:=pieces[0];
                 dUYear[m]:=StrToInt(pieces[1]);
                 dUMonth[m]:=StrToInt(pieces[2]);
                 dUDay[m]:=StrToInt(pieces[3]);
                 dUHour[m]:=StrToInt(pieces[4]);
                 dUMinute[m]:=StrToInt(pieces[5]);
                 dUSeconds[m]:=round(StrToFloat(pieces[6]));
                 dYear[m]:=StrToInt(pieces[7]);
                 dMonth[m]:=StrToInt(pieces[8]);
                 dDay[m]:=StrToInt(pieces[9]);
                 dHour[m]:=StrToInt(pieces[10]);
                 dMinute[m]:=StrToInt(pieces[11]);
                 dSeconds[m]:=round(StrToFloat(pieces[12]));
                 dCelsius[m]:=StrToFloat(pieces[13]);
                 dVolts[m]:=StrToFloat(pieces[14]);
                 dMsas[m]:=StrToFloat(pieces[15]);
                 dStatus[m]:=StrToInt(pieces[16]);
                 dMoonPhase[m]:=StrToFloat(pieces[17]);
                 dMoonElev[m]:=StrToFloat(pieces[18]);
                 dMoonIllum[m]:=StrToFloat(pieces[19]);
                 dSunElev[m]:=StrToFloat(pieces[20]);
               end;
        end;


        //Memo(Format('dMsas[m]=%0.12f dSunElev[m]=%0.12f',[dMsas[m],dSunElev[m]]));

      end else begin
        //Memo('STRLNEGTH=0 ********debug3****************');
      end;

      (* if we reach the end of the input file, proceed to write out the data of the last day before terminating *)
      //if EOF(fdata) then begin
      //  Readln(fdata, Str);
      //  writeln('EOF Str=',Str);
      //  m:=m+1;
      //  writeln(Format('EOF going to LastDay m=%d',[m]));
      //  goto LastDay;
      //end;
      //writeln('***debug4**** Checked for EOF.');

      //Slength:= lstrlen(SQM_Location); (* How long is the string in SQM_Location *)
      //Slength:= Length(SQM_Location); (* How long is the string in SQM_Location *)
      (* writeln(" string -%s- Slength = %d\n", SQM_Location, Slength); *)

      (*  Calculate the number of minutes since Local time 3PM for the time associated with this SQM record,
          implement a bug fix to eliminate a problem with daylight savings time. Use the UTC time values and correct the UTC via the
          longitude of the sample. Previously the Local Time was used, which jumped an hour at the Daylight Savings Time change
          This is the old code, commented out:
       if dHour[m]>14 then begin
        minutes_since_3pm[m]:= (dHour[m]-15)*60+ dMinute[m] + Trunc(dSeconds[m] / 60.0 + 0.5);
      end
      else begin
        minutes_since_3pm[m]:= 540+dHour[m]*60+dMinute[m]+ Trunc(dSeconds[m] / 60.+0.5);
      end;
      *)

      (* new code follows *)
      dPosNeg:=1;
      if(SQM_Long < 0.0) then
        dPosNeg:= -1;

      (* assignment to an integer will cause truncation of the remainder in the following statement, as desired *)
      dHour_Delta:= Trunc(abs(SQM_Long)/15.0 * dPosNeg);
      dShift_Hour:= dUHour[m] + dHour_Delta;

      //Memo(Format(' dPosNeg= %d',[dPosNeg]));
      //Memo(Format(' dHour_Delta= %d', [dHour_Delta]));
      //Memo(Format(' dShift_Hour= %d', [dShift_Hour]));

      if(dShift_Hour > 14) then
        minutes_since_3pm[m]:= (dShift_Hour -15) *60 + dUMinute[m] + Trunc(dUSeconds[m]/60.0+0.5)
      else
        minutes_since_3pm[m]:= 540 + dShift_Hour *60 + dUMinute[m] + Trunc(dUSeconds[m]/60.0+0.5);



      (* check whether we have reached a gap in the input data time - i.e. is this
       * data point more than the specified maximum gap length in minutes beyond the last data point? *)
      if m>0 then begin
        (* calculate the number of minutes associated with the current data point time, and compare with the previous point *)
        (* handle the special case of crossing the midnight boundary *)
        if dDay[m]=dDay[m-1] then begin
          (* if here, this new point is on the same day *)
          num_minutesA:= Round(dHour[m]*60.0 + dMinute[m]);
        end else begin
          (* if here, we have crossed the midnight boundary *)
          num_minutesA:= Round(24.*60.+dMinute[m]);
        end;

        num_minutesB:= Round(dHour[m-1]*60.+dMinute[m-1]);
        timediff:= num_minutesA-num_minutesB; (* make sure timediff is positive *)
        if timediff<0 then begin
          timediff:= timediff*-1;
        end;
        if timediff>timediff_max then begin
          (* if here, we have found a time gap in the data - consider the data so  for this day to be all that there is *)
          //writeln('Found a %d minute gap in the data just after %d-%d-%d %d:%d:%d',timediff,dYear[m-1],dMonth[m-1],dDay[m-1],dHour[m-1],dMinute[m-1],{!!!a type cast? =>} {integer(}dSeconds[m-1]);
          (* handle the case of a patch of data after a data gap during the daytime and prior to 15:00.  *)
          if dHour[m]<15 then
            (* set Start flag to 3, which we check later to loop appropriately *)
            Start:= 3;

          (* so jump into the loop which calculates the second derivatives, etc and writes out the data to the output file
           * for the data prior to this data gap *)
          goto LastDay;
        end;
      end;
      (* reset the Start flag if we are already past the first day of data and if we have gone beyond the 15 hundred hour *)
      if ((Start=2) and (dHour[m]>15)) then
          Start:= 0;
      (* reset the Start flag if we are already past the first day of data and if we have reached 15 hundred hour *)
      (* for the case of a partial day due to a data gap prior to 15:00 *)
      if ((Start=3) and (dHour[m]=15)) then
          Start:= 0;

      (* We copy data for each day into memory, beginning at 15:00 hours and going all evening, night, morning to the next afternoon *)
      (* We use a flag called Start to handle reading all the other samples acquired at the beginning of the day, samples which
       * were acquired when the hour was still 15. *)
      (* Check to see if we reached 15:00 hours on this day; we assume that the data are ordered in time sequence *)
      if ((dHour[m]=15)and(Start=0)) then begin
        (* the last sample of the previous day was m-1, so we know that the previous day has values in the arrays from 0 to m-1 *)
        LastDay:
        Last:= m-1;
        //Memo(Format('LastDay: m=%d',[m]));
        msas_Sum:= 0.0;
        msas_Count:= 0.0; (* Calculation of average Msas for the day; *)
        (* loop on the day's data *)
        //for{while} k:=0 to Pred(Last+1) { k++} do begin
        for k := 0 to Last+1-1 do begin

          (* Sun is lower than 18 degrees below the horizon and the moon is lower than 10 degrees below the horizon *)
          if ((dSunElev[k]<-18.0)and(dMoonElev[k]<-10.0)) then begin
            (* tally sum and count for msas average *)
            msas_Sum:= msas_Sum+dMsas[k];
            msas_Count:= msas_Count+1.0;
          end;
        end;
        (* we will later print out the average value for this day in those records that contributed to the average, not to all records *)
        (* That is, we will print out the average value for this day for those records when the sun was less than 18 degree below the horizon*)
        (* We assign a null value (-1.0) to all points of the sun higher than -18 degrees, and all points lacking any count values *)
        //for{while} k:=0 to Pred(Last+1) { k++} do
        for k := 0 to Last+1 do
        begin
          if ((dSunElev[k]<-18.0)and(dMoonElev[k]<-10.0)) then begin
            (* handle case of no values in the msas sum *)
            msas_Avg[k]:= -1.0;
            if msas_Count>0.0 then begin
              msas_Avg[k]:= msas_Sum / msas_Count;
            end;
          end else
            msas_Avg[k]:= -1.0;

        end;
        (* Calculate Residual Standard Error values - samples are assumed to be a constant number of minutes apart;
         * Set half_range at the program command line to specify the number of samples to consider, and given the spacing
         * between SQM measurements, the number of minutes in the sample range;
         * This fits a regression line to each point of data, and with half_range set to 9 and 5-minute sample spaceing then
         * you get a range from 45 minutes before to 45 minutes after the point, for a total of 90 minutes;
         * Program calculates the deviation from the straight line, * expressed by the sum of ((observed - expected)**2 /(expected)).
        *)
        (* initialize SS and RSE array values to zero *)
        (* use kk to track array elements of array *)

        //if ((Last+1-half_range)<0) then begin
        //  Memo(Format('Error: Stopped because of non-contiguous data set. (Last+1-half_range)=%d',[Last+1-half_range]));
        //  goto Termination;
        //end;

        for kk := 0 to 299 do begin SS[kk] := 0.0; end;
        for kk := 0 to 299 do begin  RSE[kk] := 0.0; end;

        //(* set the first half_range of RSE values a giant value ; ditto for the last half_range values *)
        for kk := 0 to half_range-1 do begin  RSE[kk] := nodata1; end;


        // first check to see if we have enough points in the current day to calculate a valid standard error statistic
        // N is the number of point in the range of the standard error calculation

        if(m < N) then begin
           Memo(Format('For date %04d-%.2d-%.2d, we only have %d data points for this day/segment and cannot calculate a valid standard error.',
             [dYear[m],dMonth[m],dDay[m],m]));
           end
        else begin

        for kk := Last+1-half_range to m-1 do begin RSE[kk] := nodata1; end;//debug orig: m-1, seems to work with m

        DOF:= Extended((half_range*2)+1-2); (* set up the degrees of freedom; we estimate two parameters, the linear regression slope and y-intercept *)
        (* loop on the sample point about which we calculate the statistic *)
        //for{while} kk:=half_range to Pred(Last+1-half_range) { kk++} do
        for kk := half_range to Last+1-half_range-1 do begin
          (* initialize sums *)
          sum_x:= 0.0;
          sum_y:= 0.0;
          sum_xy:= 0.0;
          sum_x2:= 0.0;
          sum_y2:= 0.0;
          RSE[kk]:= 0.0;
          SS[kk]:= 0.0;
          (* loop across the 2*half_range +1 values and tabulate statistics *)
          //for{while} k:=kk-half_range to Pred(kk+half_range+1) { k++} do
          for k := kk-half_range to kk+half_range+1-1 do begin
            sum_x:= sum_x+Extended(minutes_since_3pm[k]);
            sum_y:= sum_y+Extended(dMsas[k]);
            sum_xy:= sum_xy+Extended(minutes_since_3pm[k]*Extended(dMsas[k]));
            sum_x2:= sum_x2+Extended(minutes_since_3pm[k]*Extended(minutes_since_3pm[k]));
            sum_y2:= sum_y2+Extended(dMsas[k])*Extended(dMsas[k]);
            //Memo(Format(' *********************** k= %d   sum_x=%.6f, sum_y=%.6f, sum_xy=%.6f, sum_x2=%.6f, sum_y2=%.6f, dMsas[k]=%.6f',[ k, sum_x, sum_y, sum_xy, sum_x2, sum_y2, dMsas[k]]));
          end;
          mean_x:= sum_x / N;
          mean_y:= sum_y / N;
          mean_xy:= sum_xy / N;
          mean_x2:= sum_x2 / N;
          //Memo(format('mean_x=%.6f, mean_y=%.6f, mean_xy=%.6f, mean_x2=%.6f,',[mean_x, mean_y, mean_xy, mean_x2]));
          slope:= (mean_xy-(mean_x*mean_y)) / (mean_x2-(mean_x*mean_x));
          yintercept:= ((mean_x2*mean_y)-(mean_xy*mean_x)) / (mean_x2-(mean_x*mean_x));
          rcorr:= (sum_xy-sum_x*sum_y / N) / sqrt((sum_x2-(sum_x*sum_x) / N)*(sum_y2-(sum_y*sum_y) / N));
          rsqrd:= rcorr*rcorr; (* calculate means *)
          (* calculate the slope and Y-intercept of the regression line *)
          //Memo(format('kk = %d   slope=%.6f  yintercept=%.6f',[ kk, slope, yintercept]));
          (* Evaluate the regression line at all points of this interval of data to get the expected values (on the regression line) and
           * use the observed and expected values in the Residual Standard Error (RSE) calculation
          *)
          //for{while} k:=kk-half_range to Pred(kk+half_range+1) { k++} do
          for k := kk-half_range to kk+half_range+1-1 do begin
            Expected[k]:= slope*Extended(minutes_since_3pm[k])+yintercept;
            Observed[k]:= Extended(dMsas[k]);
            SS[kk]:= SS[kk]+((Observed[k]-Expected[k])*(Observed[k]-Expected[k]));
          end;
          RSE[kk]:= (Sqrt(SS[kk] / DOF))*RSE_mult; (* note that we use sqrtl here, which takes a long double argument *)
          //Memo(format('-------------------SS[kk]=%.10f, DOF=%.6f, RSE_mult=%.6f, RSE[kk]=%.10f-------------------',[SS[kk], DOF, RSE_mult, RSE[kk]]));
          (* fix up for any negative values because Expected is negative *)
          if (RSE[kk]<0.0) then
            RSE[kk]:= RSE[kk]*-1.0;

          (* fix up for any "not a number" for the case of divide by zero above *)
          if RSE[kk].IsNan then
            RSE[kk]:= nodata2;

          //Memo(Format('kk = %d  RSE=%.6f',[kk,RSE[kk]]));
        end;

        end; //End of checking (m<N)


        (* now print all this day's records to the output file *)
        for k := 0 to Last do begin
          //Memo(Format('**debug2**  for k := 0 to Last+1 do begin  k=%d Last=%d',[k, Last]));
          (* Calculate a new variable - the number of days since Jan 1, 2018 *)
          days:=get_yday(dMonth[k],dDay[k],dYear[k]); (* We actually want the number of nights since Jan 1, 2018 - that is we want to count the evening and night as part of the

           //* same "day" - actually the same "night"; So if the minutes since 15:00 hours is greater than 540 (i.e. after midnight) we
           //* subtract one day from the "days" value so those times are considered part of the previous day (i.e."night") *)
          if minutes_since_3pm[k]>=540 then
            days:= days-1;

          right_ascension:= get_right_ascension(dUYear[k],dUMonth[k],dUDay[k],dUHour[k],dUMinute[k],Round(dUSeconds[k]),SQM_Long);
          //Memo(Format(' right_ascension= %8.6f',[right_ascension]));
          (* calculate right ascension for the SQM_Location *)
          SQM_RA:= (right_ascension*15.0)*(pi / 180.0);
          SQM_Dec:= SQM_Lat*(pi / 180.0);
          Galactic_Lat:= ArcSin(sin(SQM_Dec)*sin(Dec_NGP)+cos(SQM_Dec)*cos(Dec_NGP)*cos(SQM_RA-RightAscension_NGP));
          // convert from radians to degrees
          Galactic_Lat:= Galactic_Lat*(180.0 / pi);

          //Galactic_Long:= Galactic_Long_NCP-(ArcSin((cos(SQM_Dec)*sin(SQM_RA-RightAscension_NGP) / Cos(Galactic_Lat))));
          //Galactic_Long:= Galactic_Long*(180.0 / pi); (* convert right_ascension (SQM_RA) from hours to radians *)

          YY:=  cos(SQM_Dec) * sin(SQM_RA - RightAscension_NGP);
          XX:=  (sin(SQM_Dec) * cos(Dec_NGP)) - (cos(SQM_Dec) * sin(Dec_NGP) * cos(SQM_RA - RightAscension_NGP));
          Galactic_Long:=  Galactic_Long_NCP - ArcTan2(YY,XX);
          // convert from radians to degrees
          Galactic_Long:= Galactic_Long * (180./pi);

          // Make sure that Galactic_Long is a positive number
          if(Galactic_Long < 0.0) then
              Galactic_Long:= 360. + Galactic_Long;


          (* the Declination of the SQM is its Latitude, convert it from decimal degrees to radians *)
          (* the following Equations are from Wikipedia on Celestial Coordinate Systems *)
          (* we previously set up these constants:    RightAscension_NGP,  Dec_NGP,  Galactic_Long_NCP *)
          (* convert Galactic_Lat and Galactic_Long from radians to degrees *)
          (* create Galactic_elevation angle, which we will print, along with the Galactic_Latitude *)
          //if Galactic_Lat<=0.0 then
          //  Galactic_Elevation:= 90.0+Galactic_Lat
          //else
          //  Galactic_Elevation:= 90.0-Galactic_Lat;

          J2000_days:= get_J2000(dUYear[k],dUMonth[k],dUDay[k],dUHour[k],dUMinute[k],Round(dUSeconds[k]));
          //writeln(Format('%s, %12.7f, %12.7f, %04d-%02d-%02d, %02d:%02d:%02d, %04d-%02d-%02d, %02d:%02d:%02d, %.1f, %.2f, %.2f, %1d, %.1f, %.3f, %.1f, %.3f, %04d, %f, %04d, %12.7f, %12.7f, %10.5f, %10.5f, %f,%f',[SQM_Location,SQM_Lat,SQM_Long,dUYear[k],dUMonth[k],dUDay[k],dUHour[k],dUMinute[k],Int(dUSeconds[k]),dYear[k],dMonth[k],dDay[k],dHour[k],dMinute[k],Int(dSeconds[k]),dCelsius[k],dVolts[k],dMsas[k],dStatus[k],dMoonPhase[k],dMoonElev[k],dMoonIllum[k],dSunElev[k],minutes_since_3pm[k],msas_Avg[k],days,right_ascension,Galactic_Lat,Galactic_Elevation,Galactic_Long,J2000_days,RSE[k]]));

          //Memo(Format('%s, %12.7f, %12.7f, %04d-%02d-%02d, %.2d:%.2d:%.2d, %04d-%02d-%02d, %.2d:%.2d:%.2d, %.1f, %.2f, %.2f, %1d, %.1f, %.3f, %.1f, %.3f, %04d, %f, %04d, %12.7f, %12.7f, %10.5f, %10.5f, %10.6f,%10.6f',
          //[SQM_Location,SQM_Lat,SQM_Long,dUYear[k],dUMonth[k],dUDay[k],dUHour[k],dUMinute[k],Round(dUSeconds[k]),dYear[k],dMonth[k],dDay[k],dHour[k],dMinute[k],Round(dSeconds[k]),dCelsius[k],dVolts[k],dMsas[k],dStatus[k],dMoonPhase[k],dMoonElev[k],dMoonIllum[k],dSunElev[k],minutes_since_3pm[k],msas_Avg[k],days,right_ascension,Galactic_Lat,Galactic_Elevation,Galactic_Long,J2000_days,RSE[k]]));
          (*  get the J2000 day value *)
          outstr:=LocationName+','+
                  Format('%12.7f,%12.7f,',[SQM_Lat,SQM_Long])+
                  Format('%04d-%.2d-%.2d,',[dUYear[k],dUMonth[k],dUDay[k]])+
                  Format('%.2d:%.2d:%.2d,',[dUHour[k],dUMinute[k],Round(dUSeconds[k])])+
                  Format('%04d-%.2d-%.2d,',[dYear[k],dMonth[k],dDay[k]])+
                  Format('%.2d:%.2d:%.2d,',[dHour[k],dMinute[k],Round(dSeconds[k])])+
                  Format('%.1f,',[dCelsius[k]]);
          if VoltageField>0 then
            outstr:=outstr+Format('%.2f,',[dVolts[k]])
          else
            outstr:=outstr+'-1.0,';
          outstr:=outstr+Format('%.2f,',[dMsas[k]]);
          if RecordTypeField>0 then
            outstr:=outstr+Format('%1d,',[dStatus[k]])
          else
            outstr:=outstr+'-1,';
          outstr:=outstr+
                  Format('%.1f,%.3f,%.1f,%.3f,',[dMoonPhase[k],dMoonElev[k],dMoonIllum[k],dSunElev[k]])+
                  Format('%.4d,',[minutes_since_3pm[k]])+
                  Format('%1.6f,',[msas_Avg[k]])+
                  Format('%04d,',[days])+
                  //Format('%12.7f,%12.7f,%10.5f,%10.5f,',[right_ascension,Galactic_Lat,Galactic_Elevation,Galactic_Long])+
                  Format('%12.7f,%12.7f,%10.5f,',[right_ascension,Galactic_Lat,Galactic_Long])+
                  Format('%1.6f,%.6f',[J2000_days,RSE[k]]);

          Writeln(fdataout,outstr);
          (* Note, we need to output two numbers for each of hour, minute and seconds. If only one digit is output,
             Spotfire, and other programs, will take the digit as a ten's value, insted of a one's value*)
        end;

        (* if we are at the EOF, we have already written out the last day's data, so terminate *)
        if EOF(fdata) then begin
           //Memo('ret=EOF going to Termination');
          Memo(' Reached the End of input File');
          goto Termination;
        end;

        dUYear[0]:= dUYear[m];
        dUMonth[0]:= dUMonth[m];
        dUDay[0]:= dUDay[m];
        dUHour[0]:= dUHour[m];
        dUMinute[0]:= dUMinute[m];
        dUSeconds[0]:= dUSeconds[m];
        dYear[0]:= dYear[m];
        dMonth[0]:= dMonth[m];
        dDay[0]:= dDay[m];
        dHour[0]:= dHour[m];
        dMinute[0]:= dMinute[m];
        dSeconds[0]:= dSeconds[m];
        dCelsius[0]:= dCelsius[m];
        dVolts[0]:= dVolts[m];
        dMsas[0]:= dMsas[m];
        dStatus[0]:= dStatus[m];
        dMoonPhase[0]:= dMoonPhase[m];
        dMoonElev[0]:= dMoonElev[m];
        dMoonIllum[0]:= dMoonIllum[m];
        dSunElev[0]:= dSunElev[m];
        minutes_since_3pm[0]:= minutes_since_3pm[m];
        RSE[0]:= RSE[m] / RSE_mult;
        m:= 0; (* if here, we have written out all of the day's attributes, so keep the very last record and proceed to read the next record *)
        (* m is incremented above, so set it to zero here; this avoids writing over the data we just stored at location zero *)
        if Start=3 then
          (* if here, we have a case of a partial day of data after a data gap and prior to 15:00 in the day *)
          goto ReadAnother

        else
          (* set the flag to direct all the next samples acquired in the "15" hundred hour into a new day *)
          Start:= 2;

        goto ReadAnother;
      end;

      goto ReadAnother;
      (* if here, we have reached the EOF on fscanf *)
      Termination:
      Memo(Format(' Results written to: %s',[NameOut]));
      if ErrorInputLineCounter>0 then
        Memo(Format('We found %d bad records due to missing data, and we ignored them.',[ErrorInputLineCounter]));
      Memo('Finished');
      CloseFile(fdata);
      CloseFile(fdataout);

  end;

procedure TCloudRemMilkyWay.SourceFileButtonClick(Sender: TObject);
begin

  //Clear out memo for new messages
  CloudRemMilkyWay.ProcessStatusMemo.Clear;

  SourceFileDialog.FileName:=SourceFileName;
  if SourceFileDialog.Execute then begin
    SourceFileName:=RemoveMultiSlash(SourceFileDialog.FileName);
    SourceFileEdit.Text:=SourceFileName;
  end;

  //Save directory name in registry
  vConfigurations.WriteString(Section,'SourceFileName',SourceFileName);

  FileTypeCheck;

end;

procedure TCloudRemMilkyWay.FormShow(Sender: TObject);
begin
  //For debugging/development purposes, recall last used file details.
  if ParameterCommand('-TCMR') then begin
    SourceFileName:=RemoveMultiSlash(vConfigurations.ReadString(Section, 'SourceFileName', ''));
    SourceFileEdit.Text:=SourceFileName;

    LocationName:=vConfigurations.ReadString(Section,'LocatioName','0');
    LocationEdit.Text:=LocationName;

    SQM_Lat:= StrToFloat(vConfigurations.ReadString(Section,'Latitude','0'));
    LatEdit.Text:=Format('%.7f',[SQM_Lat]);

    SQM_Long:=StrToFloat(vConfigurations.ReadString(Section,'Longitude','0'));
    LongEdit.Text:=Format('%.7f',[SQM_Long]);
  end;

  half_range:=StrToInt(vConfigurations.ReadString(Section,'half_range','9'));
  HalfRangeEdit.Text:=Format('%d',[half_range]);


  FileTypeCheck;
end;


initialization
  {$I CloudRemUnit.lrs}

end.