DPS File Concatenator - Concatenates two DPS SPS files

(* File names                                                   *)
filename="F:\\Observation Records\\AJ4CO Observatory\\DPS\\2308\\230816\\230816000000-AJ4CO DPS.sps";
filename2="F:\\Observation Records\\AJ4CO Observatory\\DPS\\2308\\230816\\230816150946-AJ4CO DPS.sps";

infile=OpenRead[filename,BinaryFormat->True];
infile2=OpenRead[filename2,BinaryFormat->True];

outfile=OpenWrite[StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>"_concatenated.sps",BinaryFormat->True];

(* Define read string function                                  *)
(* file = input file stream                                     *)
(* len = number of characters to retrieve                       *)
ReadString[file_,len_]:=
Module[{f=file,l=len},
s="";
Do[{
s=s<>BinaryRead[f, "Character8"];
},
{i,1,l}
];
s
];

(* Define proper rounding function as MM only rounds 0.5 up when the result is an even number *)
SetAttributes[Round2,Listable];
Round2[a_]:=If[FractionalPart[a]≥1/2,Ceiling[a],Floor[a]];

(* Read fixed length portion of file header                     *)
chartversion="";
chartversion=ReadString[infile,10];
chartversion2="";
chartversion2=ReadString[infile2,10];
WriteString[outfile,chartversion];

chartstart="";
chartstart=BinaryRead[infile,"Real64"];
chartstart2="";
chartstart2=BinaryRead[infile2,"Real64"];
BinaryWrite[outfile,chartstart,"Real64"];

chartend="";
chartend=BinaryRead[infile,"Real64"];
chartend2="";
chartend2=BinaryRead[infile2,"Real64"];
BinaryWrite[outfile,chartend2,"Real64"];

chartlat="";
chartlat=BinaryRead[infile,"Real64"];
chartlat2="";
chartlat2=BinaryRead[infile2,"Real64"];
BinaryWrite[outfile,chartlat,"Real64"];

chartlong="";
chartlong=BinaryRead[infile,"Real64"];
chartlong2="";
chartlong2=BinaryRead[infile2,"Real64"];
BinaryWrite[outfile,chartlong,"Real64"];

chartmax="";
chartmax=BinaryRead[infile,"Real64"];
chartmax2="";
chartmax2=BinaryRead[infile2,"Real64"];
BinaryWrite[outfile,chartmax,"Real64"];

chartmin="";
chartmin=BinaryRead[infile,"Real64"];
chartmin2="";
chartmin2=BinaryRead[infile2,"Real64"];
BinaryWrite[outfile,chartmin,"Real64"];

chartTZ="";
chartTZ=BinaryRead[infile,"Integer16"];
chartTZ2="";
chartTZ2=BinaryRead[infile2,"Integer16"];
BinaryWrite[outfile,chartTZ,"Integer16"];

chartsource="";
chartsource=ReadString[infile,10];
chartsource2="";
chartsource2=ReadString[infile2,10];
WriteString[outfile,chartsource];

chartauthor="";
chartauthor=ReadString[infile,20];
chartauthor2="";
chartauthor2=ReadString[infile2,20];
WriteString[outfile,chartauthor];

chartlocalname="";
chartlocalname=ReadString[infile,20];
chartlocalname2="";
chartlocalname2=ReadString[infile2,20];
WriteString[outfile,chartlocalname];

chartlocation="";
chartlocation=ReadString[infile,40];
chartlocation2="";
chartlocation2=ReadString[infile2,40];
WriteString[outfile,chartlocation];

chartchannels="";
chartchannels=BinaryRead[infile,"Integer16"];
chartchannels2="";
chartchannels2=BinaryRead[infile2,"Integer16"];
BinaryWrite[outfile,chartchannels,"Integer16"];

chartnotelength="";
chartnotelength=BinaryRead[infile,"Integer32"];
chartnotelength2="";
chartnotelength2=BinaryRead[infile2,"Integer32"];

chartnote="";
chartnote=ReadString[infile,chartnotelength];
chartnote2="";
chartnote2=ReadString[infile2,chartnotelength2];

(* Read number of frequency sweeps from neader notes            *)
sweepstext=StringTake[chartnote,{ StringPosition[chartnote,"SWEEPS"][[1,2]]+1 ,StringPosition[chartnote,"ÿLOWF"][[1,1]]-1}];
sweeps=ToExpression[StringTake[chartnote,{ StringPosition[chartnote,"SWEEPS"][[1,2]]+1 ,StringPosition[chartnote,"ÿLOWF"][[1,1]]-1}]];
sweeps2=ToExpression[StringTake[chartnote2,{ StringPosition[chartnote2,"SWEEPS"][[1,2]]+1 ,StringPosition[chartnote2,"ÿLOWF"][[1,1]]-1}]];

chartduration=Round2[(chartend-chartstart)*86400000]*0.001;
chartduration2=Round2[(chartend2-chartstart2)*86400000]*0.001;

sweeprate=sweeps/chartduration;
sweeprate2=sweeps2/chartduration2;

deadtime=Round2[(chartstart2-chartend)*86400000]*0.001;
deadsweeps=Round2[deadtime*sweeprate];
totalsweeps=sweeps+deadsweeps+sweeps2;

chartnote=StringReplace[chartnote,sweepstext->ToString[totalsweeps]];
chartnotelength=StringLength[chartnote];
BinaryWrite[outfile,chartnotelength,"Integer32"];
WriteString[outfile,chartnote];

Print["Chart duration = ",chartduration]
Print["Chart2 duration = ",chartduration2]
Print["Chart sweeps = ",sweeps]
Print["Chart2 sweeps = ",sweeps2]
Print["Dead time = ",deadtime]
Print["Dead sweeps = ",deadsweeps]
Print["Total sweeps = ",totalsweeps]
Print["Chart note length = ",chartnotelength]
Print["Chart note = ",chartnote]

(* Create sweep arrays, then populate and write to file         *)
(* iterating over all sweeps in one file.                       *)
indataB=Table[0,{i,1,chartchannels}];
indataA=Table[0,{i,1,chartchannels}];
Do[
{
(* Read one frequency sweep                                  *)
Do[
{
indataB[[i]]=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
indataA[[i]]=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
}
,{i,1,chartchannels}
];
(* Rotate channel 200 back down to below channel 199         *)
(* indataB=RotateLeft[indataB]; *)
(* indataA=RotateLeft[indataA]; *)
(* Write data to output file                                 *)
Do[
{
BinaryWrite[outfile,indataB[[i]],"UnsignedInteger16",ByteOrdering->+1];
BinaryWrite[outfile,indataA[[i]],"UnsignedInteger16",ByteOrdering->+1];
}
,{i,1,chartchannels}
];
(* Write EOR marker                                          *)
indataB[[1]]=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
BinaryWrite[outfile,indataB[[1]],"UnsignedInteger16",ByteOrdering->+1];
}
,{j,1,sweeps}
];

Close[infile];

Do[
{
Do[
{
BinaryWrite[outfile,0,"UnsignedInteger16",ByteOrdering->+1];
BinaryWrite[outfile,0,"UnsignedInteger16",ByteOrdering->+1];
}
,{i,1,chartchannels}
];
(* Write EOR marker                                          *)
BinaryWrite[outfile,65278,"UnsignedInteger16",ByteOrdering->+1];
}
,{j,1,deadsweeps}
];

indataB=Table[0,{i,1,chartchannels}];
indataA=Table[0,{i,1,chartchannels}];
Do[
{
(* Read one frequency sweep                                  *)
Do[
{
indataB[[i]]=BinaryRead[infile2,"UnsignedInteger16",ByteOrdering->+1];
indataA[[i]]=BinaryRead[infile2,"UnsignedInteger16",ByteOrdering->+1];
}
,{i,1,chartchannels}
];
(* Rotate channel 200 back down to below channel 199         *)
(* indataB=RotateLeft[indataB]; *)
(* indataA=RotateLeft[indataA]; *)
(* Write data to output file                                 *)
Do[
{
BinaryWrite[outfile,indataB[[i]],"UnsignedInteger16",ByteOrdering->+1];
BinaryWrite[outfile,indataA[[i]],"UnsignedInteger16",ByteOrdering->+1];
}
,{i,1,chartchannels}
];
(* Write EOR marker                                          *)
indataB[[1]]=BinaryRead[infile2,"UnsignedInteger16",ByteOrdering->+1];
BinaryWrite[outfile,indataB[[1]],"UnsignedInteger16",ByteOrdering->+1];
}
,{j,1,sweeps2}
];
Close[infile2];

Close[outfile];
Print["Done."]

DPS File Correction Array Applyer and Polarization Subtractor and Side Swapper  - De-rotates and applies a correction array to a DPS SPS file and saves as a new SPS file

(********************************************************************************************)
(*                                                                                          *)
(*  Apply correction array and write new SPS file                                           *)
(*                                                                                          *)
(*                                                                                          *)
(********************************************************************************************)

(* Open one files for input                                     *)
filename="F:\\Observation Records\\AJ4CO_DPS\\2407\\240709\\240709143000.sps";

doPolSubtract = True;

doSideSwap = False;

doPolAverage = False;


(* Define correction array to be applied *)
carrfilename="F:\\Observation Records\\AJ4CO_DPS\\DPS Correction Array 2017 01 10.csv";

outfilename=StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>" corrected using CA 2017 01 10.sps";

infile=OpenRead[filename,BinaryFormat->True];
outfile=OpenWrite[outfilename,BinaryFormat->True];
carr=Import[carrfilename];

(* Define read string function                                  *)
(* file = input file stream                                     *)
(* len = number of characters to retrieve                       *)
ReadString[file_,len_]:=
Module[{f=file,l=len},
s="";
Do[{
s=s<>BinaryRead[f, "Character8"];
},
{i,1,l}
];
s
];

(* Define proper rounding function as MM only rounds 0.5 up when the result is an even number *)
SetAttributes[Round2,Listable];
Round2[a_]:=If[FractionalPart[a]≥1/2,Ceiling[a],Floor[a]];

(* Read & wrtie fixed length portion of file header             *)
chartversion="";
chartversion=ReadString[infile,10];
WriteString[outfile,chartversion];
chartstart="";
chartstart=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartstart,"Real64"];
chartend="";
chartend=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartend,"Real64"];
chartlat="";
chartlat=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartlat,"Real64"];
chartlong="";
chartlong=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartlong,"Real64"];
chartmax="";
chartmax=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartmax,"Real64"];
chartmin="";
chartmin=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartmin,"Real64"];
chartTZ="";
chartTZ=BinaryRead[infile,"Integer16"];
BinaryWrite[outfile,chartTZ,"Integer16"];
chartsource="";
chartsource=ReadString[infile,10];
WriteString[outfile,chartsource];
chartauthor="";
chartauthor=ReadString[infile,20];
WriteString[outfile,chartauthor];
chartlocalname="";
chartlocalname=ReadString[infile,20];
WriteString[outfile,chartlocalname];
chartlocation="";
chartlocation=ReadString[infile,40];
WriteString[outfile,chartlocation];
chartchannels="";
chartchannels=BinaryRead[infile,"Integer16"];
BinaryWrite[outfile,chartchannels,"Integer16"];
chartnotelength="";
chartnotelength=BinaryRead[infile,"Integer32"];
BinaryWrite[outfile,chartnotelength,"Integer32"];
chartnote="";
chartnote=ReadString[infile,chartnotelength];
WriteString[outfile,chartnote];

(* Read number of frequency sweeps from header notes            *)
sweeps=ToExpression[StringTake[chartnote,{ StringPosition[chartnote,"SWEEPS"][[1,2]]+1 ,StringPosition[chartnote,"ÿLOWF"][[1,1]]-1}]];

(* Create sweep arrays, then populate, rotate, correct, write to file    *)
(* iterating over all sweeps in source file.                             *)
indataB=Table[0,{i,1,chartchannels}];
indataA=Table[0,{i,1,chartchannels}];
Do[
{
(* Read one frequency sweep                                  *)
Do[
{
indataB[[i]]=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
indataA[[i]]=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
}
,{i,1,chartchannels}
];
(* Rotate highest channel back down to lowest channel position         *)
indataB=RotateLeft[indataB];
indataA=RotateLeft[indataA];
(* Apply correction array                                         *)
indataB=indataB+carr[[All,2]];
indataA=indataA+carr[[All,1]];

(* Validate values and write data to output file                     *)
Do[
{
If[indataB[[i]]<0,indataB[[i]]=0];
If[indataB[[i]]>4095,indataB[[i]]=4095];
BinaryWrite[outfile,indataB[[i]],"UnsignedInteger16",ByteOrdering->+1];
If[indataA[[i]]<0,indataA[[i]]=0];
If[indataA[[i]]>4095, indataA[[i]]=4095];
BinaryWrite[outfile,indataA[[i]],"UnsignedInteger16",ByteOrdering->+1];
}
,{i,1,chartchannels}
];
(* Read & write EOR marker                                          *)
indataB[[1]]=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
BinaryWrite[outfile,indataB[[1]],"UnsignedInteger16",ByteOrdering->+1];
(* The EOR marker is two consecutive 0xFE characters, ASCII "þ" *)
(* This is seens as 65278 decimal when read as a 2-byte integer  *)
(* BinaryWrite[outfile,65278,"UnsignedInteger16",ByteOrdering->+1]; *)
}
,{j,1,sweeps}
];
Close[infile];
Close[outfile]


If[ doPolSubtract==True,
{
(********************************************************************************************)
(*                                                                                          *)
(*  Polarization Subtractor                                                                 *)
(*                                                                                          *)
(*                                                                                          *)
(********************************************************************************************)

(* Open one files for input                                     *)
filename2=outfilename;
outfilename2=StringTake[filename2,StringPosition[filename2,".sps"][[1,1]]-1]<>"_RCP-LCP_Differences.sps";
infile=OpenRead[filename2,BinaryFormat->True];
outfile=OpenWrite[outfilename2,BinaryFormat->True];

(* Define read string function                                  *)
(* file = input file stream                                     *)
(* len = number of characters to retrieve                       *)
ReadString[file_,len_]:=
Module[{f=file,l=len},
s="";
Do[{
s=s<>BinaryRead[f, "Character8"];
},
{i,1,l}
];
s
];

(* Read & wrtie fixed length portion of file header             *)
chartversion="";
chartversion=ReadString[infile,10];
WriteString[outfile,chartversion];
chartstart="";
chartstart=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartstart,"Real64"];
chartend="";
chartend=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartend,"Real64"];
chartlat="";
chartlat=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartlat,"Real64"];
chartlong="";
chartlong=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartlong,"Real64"];
chartmax="";
chartmax=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartmax,"Real64"];
chartmin="";
chartmin=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartmin,"Real64"];
chartTZ="";
chartTZ=BinaryRead[infile,"Integer16"];
BinaryWrite[outfile,chartTZ,"Integer16"];
chartsource="";
chartsource=ReadString[infile,10];
WriteString[outfile,chartsource];
chartauthor="";
chartauthor=ReadString[infile,20];
WriteString[outfile,chartauthor];
chartlocalname="";
chartlocalname=ReadString[infile,20];
WriteString[outfile,chartlocalname];
chartlocation="";
chartlocation=ReadString[infile,40];
WriteString[outfile,chartlocation];
chartchannels="";
chartchannels=BinaryRead[infile,"Integer16"];
BinaryWrite[outfile,chartchannels,"Integer16"];
chartnotelength="";
chartnotelength=BinaryRead[infile,"Integer32"];
BinaryWrite[outfile,chartnotelength,"Integer32"];
chartnote="";
chartnote=ReadString[infile,chartnotelength];
WriteString[outfile,chartnote];

(* Read number of frequency sweeps from neader notes            *)
sweeps=ToExpression[StringTake[chartnote,{ StringPosition[chartnote,"SWEEPS"][[1,2]]+1 ,StringPosition[chartnote,"ÿLOWF"][[1,1]]-1}]];

(* Read data, write to file in swapped order                    *)
Do[
{
Do[
{
(* Read Side B data                                       *)
indataB=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
(* Read Side A data                                       *)
indataA=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
A=indataA;
B=indataB;
(* Subtract LCP value from RCP value and store in RCP position *)
(* Multiply difference by 5 to make it easier to see in the spectrogram display *)
(* Keep values in 12-bit range of 0 to 4095 inclusive          *)
If[A-B≥0,If[5*(A-B)>4095,indataA=4095,indataA=5*(A-B)],indataA=0];
(* Subtract RCP value from LCP value and store in LCP position *)
If[B-A≥0,If[5*(B-A)>4095,indataB=4095,indataB=5*(B-A)],indataB=0];
(* Write Side B data                                      *)
BinaryWrite[outfile,indataB,"UnsignedInteger16",ByteOrdering->+1];
(* Write Side A data                                      *)
BinaryWrite[outfile,indataA,"UnsignedInteger16",ByteOrdering->+1];
}
,{i,1,chartchannels}
];
(* Read & write EOR marker                                  *)
indataB=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
BinaryWrite[outfile,indataB,"UnsignedInteger16",ByteOrdering->+1];
}
,{j,1,sweeps}
];
Close[infile];
Close[outfile]

},
{}
]


If[ doPolAverage==True,
{
(********************************************************************************************)
(*                                                                                          *)
(*  Polarization Averager                                                                   *)
(*                                                                                          *)
(*                                                                                          *)
(********************************************************************************************)

(* Open one files for input                                     *)
filename4=outfilename;
outfilename4=StringTake[filename4,StringPosition[filename4,".sps"][[1,1]]-1]<>"_RCP-LCP_Averages.sps";
infile=OpenRead[filename4,BinaryFormat->True];
outfile=OpenWrite[outfilename4,BinaryFormat->True];

(* Define read string function                                  *)
(* file = input file stream                                     *)
(* len = number of characters to retrieve                       *)
ReadString[file_,len_]:=
Module[{f=file,l=len},
s="";
Do[{
s=s<>BinaryRead[f, "Character8"];
},
{i,1,l}
];
s
];

(* Read & wrtie fixed length portion of file header             *)
chartversion="";
chartversion=ReadString[infile,10];
WriteString[outfile,chartversion];
chartstart="";
chartstart=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartstart,"Real64"];
chartend="";
chartend=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartend,"Real64"];
chartlat="";
chartlat=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartlat,"Real64"];
chartlong="";
chartlong=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartlong,"Real64"];
chartmax="";
chartmax=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartmax,"Real64"];
chartmin="";
chartmin=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartmin,"Real64"];
chartTZ="";
chartTZ=BinaryRead[infile,"Integer16"];
BinaryWrite[outfile,chartTZ,"Integer16"];
chartsource="";
chartsource=ReadString[infile,10];
WriteString[outfile,chartsource];
chartauthor="";
chartauthor=ReadString[infile,20];
WriteString[outfile,chartauthor];
chartlocalname="";
chartlocalname=ReadString[infile,20];
WriteString[outfile,chartlocalname];
chartlocation="";
chartlocation=ReadString[infile,40];
WriteString[outfile,chartlocation];
chartchannels="";
chartchannels=BinaryRead[infile,"Integer16"];
BinaryWrite[outfile,chartchannels,"Integer16"];
chartnotelength="";
chartnotelength=BinaryRead[infile,"Integer32"];
BinaryWrite[outfile,chartnotelength,"Integer32"];
chartnote="";
chartnote=ReadString[infile,chartnotelength];
WriteString[outfile,chartnote];

(* Read number of frequency sweeps from neader notes            *)
sweeps=ToExpression[StringTake[chartnote,{ StringPosition[chartnote,"SWEEPS"][[1,2]]+1 ,StringPosition[chartnote,"ÿLOWF"][[1,1]]-1}]];

(* Read data, write to file in swapped order                    *)
Do[
{
Do[
{
(* Read Side B data                                       *)
indataB=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
(* Read Side A data                                       *)
indataA=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
A=indataA;
B=indataB;
(* Average RCP and LCP          *)
indataA=Round2[(A+B)/2];
indataB=indataA;
(* Write Side B data                                      *)
BinaryWrite[outfile,indataB,"UnsignedInteger16",ByteOrdering->+1];
(* Write Side A data                                      *)
BinaryWrite[outfile,indataA,"UnsignedInteger16",ByteOrdering->+1];
}
,{i,1,chartchannels}
];
(* Read & write EOR marker                                  *)
indataB=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
BinaryWrite[outfile,indataB,"UnsignedInteger16",ByteOrdering->+1];
}
,{j,1,sweeps}
];
Close[infile];
Close[outfile]

},
{}
]


If[doSideSwap==True,
{
(********************************************************************************************)
(*                                                                                          *)
(*  Swap Side A with Side B data and write new SPS file                                     *)
(*                                                                                          *)
(*                                                                                          *)
(********************************************************************************************)

(* Open one files for input                                     *)
filename3=outfilename;
outfilename3=StringTake[filename3,StringPosition[filename3,".sps"][[1,1]]-1]<>"_rev_pol.sps";
infile=OpenRead[filename3,BinaryFormat->True];
outfile=OpenWrite[outfilename3,BinaryFormat->True];

(* Define read string function                                  *)
(* file = input file stream                                     *)
(* len = number of characters to retrieve                       *)
ReadString[file_,len_]:=
Module[{f=file,l=len},
s="";
Do[{
s=s<>BinaryRead[f, "Character8"];
},
{i,1,l}
];
s
];

(* Read & wrtie fixed length portion of file header             *)
chartversion="";
chartversion=ReadString[infile,10];
WriteString[outfile,chartversion];
chartstart="";
chartstart=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartstart,"Real64"];
chartend="";
chartend=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartend,"Real64"];
chartlat="";
chartlat=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartlat,"Real64"];
chartlong="";
chartlong=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartlong,"Real64"];
chartmax="";
chartmax=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartmax,"Real64"];
chartmin="";
chartmin=BinaryRead[infile,"Real64"];
BinaryWrite[outfile,chartmin,"Real64"];
chartTZ="";
chartTZ=BinaryRead[infile,"Integer16"];
BinaryWrite[outfile,chartTZ,"Integer16"];
chartsource="";
chartsource=ReadString[infile,10];
WriteString[outfile,chartsource];
chartauthor="";
chartauthor=ReadString[infile,20];
WriteString[outfile,chartauthor];
chartlocalname="";
chartlocalname=ReadString[infile,20];
WriteString[outfile,chartlocalname];
chartlocation="";
chartlocation=ReadString[infile,40];
WriteString[outfile,chartlocation];
chartchannels="";
chartchannels=BinaryRead[infile,"Integer16"];
BinaryWrite[outfile,chartchannels,"Integer16"];
chartnotelength="";
chartnotelength=BinaryRead[infile,"Integer32"];
BinaryWrite[outfile,chartnotelength,"Integer32"];
chartnote="";
chartnote=ReadString[infile,chartnotelength];
WriteString[outfile,chartnote];

(* Read number of frequency sweeps from neader notes            *)
sweeps=ToExpression[StringTake[chartnote,{ StringPosition[chartnote,"SWEEPS"][[1,2]]+1 ,StringPosition[chartnote,"ÿLOWF"][[1,1]]-1}]];

(* Read data, write to file in swapped order                    *)
Do[
{
Do[
{
(* Read Side B data                                       *)
indataB=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
(* Read Side A data                                       *)
indataA=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
(* Write Side A data in Side B position                   *)
BinaryWrite[outfile,indataA,"UnsignedInteger16",ByteOrdering->+1];
(* Write Side B data in Side A position                   *)
BinaryWrite[outfile,indataB,"UnsignedInteger16",ByteOrdering->+1];
}
,{i,1,chartchannels}
];
(* Read & write EOR marker                                  *)
indataB=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
BinaryWrite[outfile,indataB,"UnsignedInteger16",ByteOrdering->+1];
}
,{j,1,sweeps}
];
Close[infile];
Close[outfile]
Print["Done."]
},
{}
]

DPS Data Analyzer - Multiple Files - Full analysis of DPS data with correction arrays

(********************************************************************************************)
(*                                                                                          *)
(*  Takes an SPS file and plots Side A and Side B spectrograms and average continua.        *)
(*  Saves a PDF report in the same folder as the source data.                               *)
(*  Saves PNG files of the plots in the same folder as the source data.                     *)
(*                                                                                          *)
(********************************************************************************************)

(* Specify the RSP data root directory *)
dataroot="T:\\AJ4CO Archive\\AJ4CO Observatory\\DPS\\1605\\160512";

(* Specify output directory for all image files *) 
outputdir="T:\\AJ4CO Archive\\AJ4CO Observatory\\DPS\\Individual Spectrograms\\2016 05";

(* Specify spectrgram plot width in seconds *)
plotduration = 86400;

(* Set file name of correction array                            *)
cafilename="R:\\Observation Records\\AJ4CO_DPS\\DPS Correction Array 2017 01 10.csv";

(* Set color offset and color gain for raw RCP and LCP spectrograms *)
rawcoloroffset=1975;
rawcolorgain=1.95;

(* Set color offset and color gain for corrected RCP and LCP spectrograms *)
corrcoloroffset=2075;
corrcolorgain=5.0;

(* Set color offset and color gain for stoke -V and +V spectrograms *)
stokescoloroffset=50;
stokescolorgain=2.0;

(* Set output chart width in pixels *)
chartwidth=1440;

(* Set color file path and name for RCP & LCP spectrograms      *)
cpcolorfile="N:\\Documents\\0 RA\\Software\\Radio-Sky Spectrograph\\Color File Maker\\AJ4CO-Rainbow.txt";

(* Set number of bits in data                                   *)
nbits=12;

(* Define proper rounding function as MM only rounds 0.5 up when the result is an even number *)
SetAttributes[Round2,Listable];
Round2[a_]:=If[FractionalPart[a]≥1/2,Ceiling[a],Floor[a]];

(* Define read string function                                  *)
(* file = input file stream                                     *)
(* len = number of characters to retrieve                       *)
ReadString[file_,len_]:=
Module[{f=file,l=len},
s="";
Do[{
s=s<>BinaryRead[f, "Character8"];
},
{i,1,l}
];
s
];

(* Print start time stamp *)
Print["Processing started on "<>DateString[]]

(* Find all SPS files within the root directory's structure *)
fileext=".sps";
SetDirectory[dataroot];
Files=FileNames["*"<>fileext,{"*"},Infinity];
(* count the number of files found *)
filecount=Dimensions[Files][[1]];
(* log the number of files found *)
Print["Found "<>ToString[filecount]<>" "<>fileext<>" files to process."]
Print[]

(* Open and read correction array                               *)
corrarray=Transpose[Import[cafilename,"Data"]];

(* Open and read color file data for circular pol spectrograms. *)
(* This produces array "ColorsInT" used by ArrayPlot[]'s color  *)
(* function specification.                                      *)
ColorInFile=Import[cpcolorfile,"Data"];
(* first element of color data                                  *)
colorlo=1;
(* number of color data levels, usually 4096                    *)
colorhi=Dimensions[ColorInFile][[1]];
(* Extract 24-bit RGB color triplets from color file data       *)
(* Color file data has 4096 levels, each level has format:      *)
(* Color file data value = Blue*65536 + Green*256 + Red         *)
(* where R, G, and B are 8-bit integers                         *)
ColorsIn ={
(* red component *)
FractionalPart[FractionalPart[ColorInFile[[All,1]]/65536]*65536/256]*256,
(* green component *)
IntegerPart[FractionalPart[ColorInFile[[All,1]]/65536]*65536/256],
(* blue component *)
IntegerPart[ColorInFile[[All,1]]/65536]
};
(* Transpose color triplet array to match RCP & LCP arrays      *)
ColorsInT=Transpose[ColorsIn];
(* Scale color triplets to Mathematica's color domain of 0 to 1 *)
ColorsInT=ColorsInT/255;


(* Cycle over all found files  *)
Do [
{

(* Create output filenames *)
trimmedfilename=StringReplace[StringDrop[StringDrop[Files[[z]],StringPosition[Files[[z]],"\\"][[Length[StringPosition[Files[[z]],"\\"]],1]]],-4]," "->"_"];
spgmAfilename=outputdir<>"\\"<>trimmedfilename<>"_24_Hour_Spectrogram_RCP_Raw.png";
spgmBfilename=outputdir<>"\\"<>trimmedfilename<>"_24_Hour_Spectrogram_LCP_Raw.png";
spgmAfilenamecorr=outputdir<>"\\"<>trimmedfilename<>"_24_Hour_Spectrogram_RCP_Corrected.png";
spgmBfilenamecorr=outputdir<>"\\"<>trimmedfilename<>"_24_Hour_Spectrogram_LCP_Corrected.png";
spgmAfilenamecorrstokes=outputdir<>"\\"<>trimmedfilename<>"_24_Hour_Spectrogram_RCP_Corrected_Stokes.png";
spgmBfilenamecorrstokes=outputdir<>"\\"<>trimmedfilename<>"_24_Hour_Spectrogram_LCP_Corrected_Stokes.png";

(* Open one SPS file for input                                  *)
infile=OpenRead[Files[[z]],BinaryFormat->True];
Print[z,"  Processing "<>Files[[z]]<> " on "<>DateString[]];
Print[];
(* Read fixed length portion of file header                     *)
chartversion="";
chartversion=ReadString[infile,10];
chartstart="";
chartstart=BinaryRead[infile,"Real64"];
chartend="";
chartend=BinaryRead[infile,"Real64"];
chartlat="";
chartlat=BinaryRead[infile,"Real64"];
chartlong="";
chartlong=BinaryRead[infile,"Real64"];
chartmax="";
chartmax=BinaryRead[infile,"Real64"];
chartmin="";
chartmin=BinaryRead[infile,"Real64"];
chartTZ="";
chartTZ=BinaryRead[infile,"Integer16"];
chartsource="";
chartsource=ReadString[infile,10];
chartauthor="";
chartauthor=ReadString[infile,20];
chartlocalname="";
chartlocalname=ReadString[infile,20];
chartlocation="";
chartlocation=ReadString[infile,40];
chartchannels="";
chartchannels=BinaryRead[infile,"Integer16"];
chartnotelength="";
chartnotelength=BinaryRead[infile,"Integer32"];
chartnote="";
chartnote=ReadString[infile,chartnotelength];

(* Calculate parameters stored in header info and format data   *)
(* for display.                                                 *)
chartduration=Round2[(chartend-chartstart)*86400000]*0.001;
chartstartdate=DateList[(chartstart-2)*86400];
chartenddate=DateList[(chartend-2)*86400];
chartstartF=DateString[chartstartdate, "Day"]<>"-"<>DateString[chartstartdate, "MonthName"]<>"-"<>DateString[chartstartdate, "Year"]<>"  "<>DateString[chartstartdate, "Hour24"]<>":"<>DateString[chartstartdate, "Minute"]<>":"<>DateString[chartstartdate, "Second"]<>"."<>DateString[chartstartdate, "Millisecond"]<>" UTC";
chartendF=DateString[chartenddate, "Day"]<>"-"<>DateString[chartenddate, "MonthName"]<>"-"<>DateString[chartenddate, "Year"]<>"  "<>DateString[chartenddate, "Hour24"]<>":"<>DateString[chartenddate, "Minute"]<>":"<>DateString[chartenddate, "Second"]<>"."<>DateString[chartenddate, "Millisecond"]<>" UTC";
chartversionF=ToString[FromDigits[StringTake[chartversion,5]]]<>"."<>ToString[FromDigits[StringTake[chartversion,{6,7}]]]<>"."<>ToString[FromDigits[StringTake[chartversion,{8,10}]]];
freqlo=ToExpression[StringTake[chartnote,{ StringPosition[chartnote,"LOWF"][[1,2]]+1 ,StringPosition[chartnote,"LOWF"][[1,2]]+8}]]*10^-6;
freqhi=ToExpression[StringTake[chartnote,{ StringPosition[chartnote,"HIF"][[1,2]]+1 ,StringPosition[chartnote,"HIF"][[1,2]]+8}]]*10^-6;
channelsep=N[(freqhi-freqlo)/(chartchannels+1)*10^3]; (* channels+1 to account for RSS channel allocation bug *)
freqhi=N[freqhi-channelsep*10^-3]; (* account for RSS channel allocation bug *)
channeldwelltime=Round2[chartduration/(sweeps*chartchannels)*10^6];
sweeps=ToExpression[StringTake[chartnote,{ StringPosition[chartnote,"SWEEPS"][[1,2]]+1 ,StringPosition[chartnote,"ÿLOWF"][[1,1]]-1}]];
chartsweeprate=sweeps/chartduration;
If[StringPosition[chartsource,FromCharacterCode[0]][[1,1]]==1,chartsource="< not specified >"];
chartsweeps=Round2[chartsweeprate*plotduration];
shrinkratio=chartsweeps/chartwidth;
startoffset=Floor[86400*(chartstart-IntegerPart[chartstart])*chartsweeprate];

(* Read SPS data file into dataAraw and dataBraw matrices       *)
(* Iterate read function over sweeps and over the number of     *)
(* chart channels * 2, plus one EOR marker per sweep.           *)
(* NOTE: SPS data file ADC values are big endian (the byte      *)
(* order is reversed; i.e., MSB is the low byte)                *)
dataAraw=Table[0,{i,1,chartsweeps},{j,1,chartchannels}];
dataBraw=Table[0,{i,1,chartsweeps},{j,1,chartchannels}];
Do[
{
(* Read one input data record                                 *)
Do[
{
(* Read Side B data (side B comes first in the SPS data)    *)
dataBraw[[i,j]]=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
(* Read Side A data                                         *)
dataAraw[[i,j]]=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
}
(* cycle over 2 x number of channels *)
,{j,1,chartchannels}
];
(* Read EOR marker                                            *)
dummy=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
}
(* cycle over number of input sweeps to be averaged            *)
,{i,1+startoffset,sweeps+startoffset} 
];
Close[infile];


(*  RSS records channel 300 rolled over above channel 1.              *)
(*  These two lines shift channel 300 data back to where it belongs.  *)
dataAraw=RotateLeft[dataAraw,{0,1}];
dataBraw=RotateLeft[dataBraw,{0,1}];


(* Apply correction arrays to raw data                           *)
dataAcorr=dataAraw;
dataBcorr=dataBraw;
Do[dataAcorr[[i]]=dataAraw[[i]]+corrarray[[1]],{i,1,chartsweeps}];
Do[dataBcorr[[i]]=dataBraw[[i]]+corrarray[[2]],{i,1,chartsweeps}];


(* Do Stokes parameter subtractions on crrected data                    *)
dataAstokes=dataAcorr;
dataBstokes=dataBcorr;
Do[
{
A=dataAcorr[[i,j]];
B=dataBcorr[[i,j]];
Ax=0;
Bx=0;
(* Subtract LCP value from RCP value and store in RCP position *)
(* Multiply difference by 5 to make it easier to see in the spectrogram display *)
(* Keep values in 12-bit range of 0 to 4095 inclusive          *)
If[A-B≥0,If[5*(A-B)>4095,Ax=4095,Ax=5*(A-B)],Ax=0];
(* Subtract RCP value from LCP value and store in LCP position *)
If[B-A≥0,If[5*(B-A)>4095,Bx=4095,Bx=5*(B-A)],Bx=0];
(* Write Side A data                                      *)
dataAstokes[[i,j]]=Ax;
(* Write Side B data                                      *)
dataBstokes[[i,j]]=Bx;
}
,{i,1,chartsweeps},{j,1,chartchannels}];


(* Decimate raw data to fit within the spectrogram image width   *)
dataArawdecimated=Table[0,{i,1,chartwidth},{j,1,chartchannels}];
dataBrawdecimated=Table[0,{i,1,chartwidth},{j,1,chartchannels}];
Do[{
dataArawdecimated[[i,j]]=Mean[dataAraw[[Round2[1+(i-1)*shrinkratio];;Round2[i*shrinkratio],j]]];
dataBrawdecimated[[i,j]]=Mean[dataBraw[[Round2[1+(i-1)*shrinkratio];;Round2[i*shrinkratio],j]]];
},{j,1, chartchannels},{i,1, chartwidth}];


(* Decimate corrected data to fit within the spectrogram image width   *)
dataAcorrdecimated=Table[0,{i,1,chartwidth},{j,1,chartchannels}];
dataBcorrdecimated=Table[0,{i,1,chartwidth},{j,1,chartchannels}];
Do[{
dataAcorrdecimated[[i,j]]=Mean[dataAcorr[[Round2[1+(i-1)*shrinkratio];;Round2[i*shrinkratio],j]]];
dataBcorrdecimated[[i,j]]=Mean[dataBcorr[[Round2[1+(i-1)*shrinkratio];;Round2[i*shrinkratio],j]]];
},{j,1, chartchannels},{i,1, chartwidth}];


(* Decimate stokes data to fit within the spectrogram image width   *)
dataAstokesdecimated=Table[0,{i,1,chartwidth},{j,1,chartchannels}];
dataBstokesdecimated=Table[0,{i,1,chartwidth},{j,1,chartchannels}];
Do[{
dataAstokesdecimated[[i,j]]=Mean[dataAstokes[[Round2[1+(i-1)*shrinkratio];;Round2[i*shrinkratio],j]]];
dataBstokesdecimated[[i,j]]=Mean[dataBstokes[[Round2[1+(i-1)*shrinkratio];;Round2[i*shrinkratio],j]]];
},{j,1, chartchannels},{i,1, chartwidth}];


(* Apply color offset and color gain to raw data, just like RSS *) 
dataArawdecimatedadj=Round2[(dataArawdecimated-rawcoloroffset)*rawcolorgain];
dataBrawdecimatedadj=Round2[(dataBrawdecimated-rawcoloroffset)*rawcolorgain];
(* Coerce color-adjusted data to fit wihtin the n-bit range     *)
datalo=1;
datahi=2^nbits-1;  
(* Coerce Side A data                                          *)
Do[{
If[
dataArawdecimatedadj[[i,j]]≤datahi,
{If[
dataArawdecimatedadj[[i,j]]≥datalo,
{dataArawdecimatedadj[[i,j]]=dataArawdecimatedadj[[i,j]]},
{dataArawdecimatedadj[[i,j]]=datalo}
]},{dataArawdecimatedadj[[i,j]]=datahi}];
},{i,1,chartwidth},{j,1,chartchannels}];
(* Coerce Side B data                                          *)
Do[{
If[
dataBrawdecimatedadj[[i,j]]≤datahi,
{If[
dataBrawdecimatedadj[[i,j]]≥datalo,
{dataBrawdecimatedadj[[i,j]]=dataBrawdecimatedadj[[i,j]]},
{dataBrawdecimatedadj[[i,j]]=datalo}
]},{dataBrawdecimatedadj[[i,j]]=datahi}]
},{i,1,chartwidth},{j,1,chartchannels}];


(* Apply color offset and color gain to corrected data, just like RSS *) 
dataAcorrdecimatedadj=Round2[(dataAcorrdecimated-corrcoloroffset)*corrcolorgain];
dataBcorrdecimatedadj=Round2[(dataBcorrdecimated-corrcoloroffset)*corrcolorgain];
(* Coerce color-adjusted data to fit wihtin the n-bit range     *)
datalo=1;
datahi=2^nbits-1;  
(* Coerce Side A data                                          *)
Do[{
If[
dataAcorrdecimatedadj[[i,j]]≤datahi,
{If[
dataAcorrdecimatedadj[[i,j]]≥datalo,
{dataAcorrdecimatedadj[[i,j]]=dataAcorrdecimatedadj[[i,j]]},
{dataAcorrdecimatedadj[[i,j]]=datalo}
]},{dataAcorrdecimatedadj[[i,j]]=datahi}];
},{i,1,chartwidth},{j,1,chartchannels}];
(* Coerce Side B data                                          *)
Do[{
If[
dataBcorrdecimatedadj[[i,j]]≤datahi,
{If[
dataBcorrdecimatedadj[[i,j]]≥datalo,
{dataBcorrdecimatedadj[[i,j]]=dataBcorrdecimatedadj[[i,j]]},
{dataBcorrdecimatedadj[[i,j]]=datalo}
]},{dataBcorrdecimatedadj[[i,j]]=datahi}]
},{i,1,chartwidth},{j,1,chartchannels}];


(* Apply color offset and color gain to stokes data, just like RSS *) 
dataAstokesdecimatedadj=Round2[(dataAstokesdecimated-stokescoloroffset)*stokescolorgain];
dataBstokesdecimatedadj=Round2[(dataBstokesdecimated-stokescoloroffset)*stokescolorgain];
(* Coerce color-adjusted data to fit wihtin the n-bit range     *)
datalo=1;
datahi=2^nbits-1;  
(* Coerce Side A data                                          *)
Do[{
If[
dataAstokesdecimatedadj[[i,j]]≤datahi,
{If[
dataAstokesdecimatedadj[[i,j]]≥datalo,
{dataAstokesdecimatedadj[[i,j]]=dataAstokesdecimatedadj[[i,j]]},
{dataAstokesdecimatedadj[[i,j]]=datalo}
]},{dataAstokesdecimatedadj[[i,j]]=datahi}];
},{i,1,chartwidth},{j,1,chartchannels}];
(* Coerce Side B data                                          *)
Do[{
If[
dataBstokesdecimatedadj[[i,j]]≤datahi,
{If[
dataBstokesdecimatedadj[[i,j]]≥datalo,
{dataBstokesdecimatedadj[[i,j]]=dataBstokesdecimatedadj[[i,j]]},
{dataBstokesdecimatedadj[[i,j]]=datalo}
]},{dataBstokesdecimatedadj[[i,j]]=datahi}]
},{i,1,chartwidth},{j,1,chartchannels}];


plot1=ArrayPlot[Transpose[dataArawdecimatedadj],
PixelConstrained->True,
ImageSize->{chartwidth,chartchannels},
ColorFunctionScaling->False,
ColorFunction->Function[{z},RGBColor[ColorsInT[[z]]]],
Background->None,
Frame->False
];

plot2=ArrayPlot[Transpose[dataBrawdecimatedadj],
PixelConstrained->True,
ImageSize->{chartwidth,chartchannels},
ColorFunctionScaling->False,
ColorFunction->Function[{z},RGBColor[ColorsInT[[z]]]],
Background->None,
Frame->False
];

plot3=ArrayPlot[Transpose[dataAcorrdecimatedadj],
PixelConstrained->True,
ImageSize->{chartwidth,chartchannels},
ColorFunctionScaling->False,
ColorFunction->Function[{z},RGBColor[ColorsInT[[z]]]],
Background->None,
Frame->False
];

plot4=ArrayPlot[Transpose[dataBcorrdecimatedadj],
PixelConstrained->True,
ImageSize->{chartwidth,chartchannels},
ColorFunctionScaling->False,
ColorFunction->Function[{z},RGBColor[ColorsInT[[z]]]],
Background->None,
Frame->False
];

plot5=ArrayPlot[Transpose[dataAstokesdecimatedadj],
PixelConstrained->True,
ImageSize->{chartwidth,chartchannels},
ColorFunctionScaling->False,
ColorFunction->Function[{z},RGBColor[ColorsInT[[z]]]],
Background->None,
Frame->False
];

plot6=ArrayPlot[Transpose[dataBstokesdecimatedadj],
PixelConstrained->True,
ImageSize->{chartwidth,chartchannels},
ColorFunctionScaling->False,
ColorFunction->Function[{z},RGBColor[ColorsInT[[z]]]],
Background->None,
Frame->False
];

(* export PNG files *)
Export[spgmAfilename,plot1];
Export[spgmBfilename,plot2];
Export[spgmAfilenamecorr,plot3];
Export[spgmBfilenamecorr,plot4];
Export[spgmAfilenamecorrstokes,plot5];
Export[spgmBfilenamecorrstokes,plot6];


},{z,1,filecount}
];

(* Print ens time stamp *)
Print["Processing Complete on "<>DateString[]]

chartstart

chartstart*86400

chartend*86400

(chartend-chartstart)*86400

86400*(chartstart-IntegerPart[chartstart])

86400*(chartend-IntegerPart[chartend])

Files

Floor[86400*(chartstart-IntegerPart[chartstart])*chartsweeprate]

Floor

DPS Data Analyzer - Multiple Files - Full analysis of DPS data with correction arrays (newer version of MM)

(********************************************************************************************)
(*                                                                                          *)
(*  Takes an SPS file and plots Side A and Side B spectrograms and average continua.       *)
(*  Saves a PDF report in the same folder as the source data.                               *)
(*  Saves PNG files of the plots in the same folder as the source data.                     *)
(*                                                                                          *)
(********************************************************************************************)

(* Specify the RSP data root directory *)
dataroot="F:\\Observation Records\\AJ4CO Observatory\\DPS\\2109";

(* Specify output directory for all image files *) 
outputdir=dataroot;

(* Specify spectrgram plot width in seconds *)
plotduration = 86400;

(* Set file name of correction array                            *)
cafilename="F:\\Observation Records\\AJ4CO_DPS\\DPS Correction Array 2017 01 10.csv";

(* Set color offset and color gain for raw RCP and LCP spectrograms *)
rawcoloroffset=1975;
rawcolorgain=1.95;

(* Set color offset and color gain for corrected RCP and LCP spectrograms *)
corrcoloroffset=2100;
corrcolorgain=5.0;

(* Set color offset and color gain for stoke -V and +V spectrograms *)
stokescoloroffset=50;
stokescolorgain=2.0;

(* Set output chart width in pixels *)
chartwidth=1440;

(* Set color file path and name for RCP & LCP spectrograms      *)
cpcolorfile="\\\\ra15\\webroot\\AJ4CO Software Archive\\Color File Maker\\AJ4CO-Rainbow.txt";

(* Set number of bits in data                                   *)
nbits=12;

(* Define proper rounding function as MM only rounds 0.5 up when the result is an even number *)
SetAttributes[Round2,Listable];
Round2[a_]:=If[FractionalPart[a]≥1/2,Ceiling[a],Floor[a]];

(* Define read string function                                  *)
(* file = input file stream                                     *)
(* len = number of characters to retrieve                       *)
ReadMyString[file_,len_]:=
Module[{f=file,l=len},
s="";
Do[{
s=s<>BinaryRead[f, "Character8"];
},
{i,1,l}
];
s
];

(* Print start time stamp *)
Print["Processing started on "<>DateString[]]

(* Find all SPS files within the root directory's structure *)
fileext=".sps";
SetDirectory[dataroot];
Files=FileNames["*"<>fileext,{"*"},Infinity];
(* count the number of files found *)
filecount=Dimensions[Files][[1]];
(* log the number of files found *)
Print["Found "<>ToString[filecount]<>" "<>fileext<>" files to process."]
Print[]

(* Open and read correction array                               *)
corrarray=Transpose[Import[cafilename,"Data"]];

(* Open and read color file data for circular pol spectrograms. *)
(* This produces array "ColorsInT" used by ArrayPlot[]'s color  *)
(* function specification.                                      *)
ColorInFile=Import[cpcolorfile,"Data"];
(* first element of color data                                  *)
colorlo=1;
(* number of color data levels, usually 4096                    *)
colorhi=Dimensions[ColorInFile][[1]];
(* Extract 24-bit RGB color triplets from color file data       *)
(* Color file data has 4096 levels, each level has format:      *)
(* Color file data value = Blue*65536 + Green*256 + Red         *)
(* where R, G, and B are 8-bit integers                         *)
ColorsIn ={
(* red component *)
FractionalPart[FractionalPart[ColorInFile[[All,1]]/65536]*65536/256]*256,
(* green component *)
IntegerPart[FractionalPart[ColorInFile[[All,1]]/65536]*65536/256],
(* blue component *)
IntegerPart[ColorInFile[[All,1]]/65536]
};
(* Transpose color triplet array to match RCP & LCP arrays      *)
ColorsInT=Transpose[ColorsIn];
(* Scale color triplets to Mathematica's color domain of 0 to 1 *)
ColorsInT=ColorsInT/255;


(* Cycle over all found files  *)
Do [
{

(* Create output filenames *)
trimmedfilename=StringReplace[StringDrop[StringDrop[Files[[z]],StringPosition[Files[[z]],"\\"][[Length[StringPosition[Files[[z]],"\\"]],1]]],-4]," "->"_"];
spgmAfilename=outputdir<>"\\"<>trimmedfilename<>"_24_Hour_Spectrogram_RCP_Raw.png";
spgmBfilename=outputdir<>"\\"<>trimmedfilename<>"_24_Hour_Spectrogram_LCP_Raw.png";
spgmAfilenamecorr=outputdir<>"\\"<>trimmedfilename<>"_24_Hour_Spectrogram_RCP_Corrected.png";
spgmBfilenamecorr=outputdir<>"\\"<>trimmedfilename<>"_24_Hour_Spectrogram_LCP_Corrected.png";
spgmAfilenamecorrstokes=outputdir<>"\\"<>trimmedfilename<>"_24_Hour_Spectrogram_RCP_Corrected_Stokes.png";
spgmBfilenamecorrstokes=outputdir<>"\\"<>trimmedfilename<>"_24_Hour_Spectrogram_LCP_Corrected_Stokes.png";

(* If all six png files have already been written, skip to next iteration *)
If[FileExistsQ[spgmBfilenamecorrstokes]==True,
{Print[Files[[z]]," already processed, skipping."];
Continue[];
}];

(* Open one SPS file for input                                  *)
infile=OpenRead[Files[[z]],BinaryFormat->True];
Print[z,"  Processing "<>Files[[z]]<> " on "<>DateString[]];
Print[];
(* Read fixed length portion of file header                     *)
chartversion="";
chartversion=ReadMyString[infile,10];
chartstart="";
chartstart=BinaryRead[infile,"Real64"];
chartend="";
chartend=BinaryRead[infile,"Real64"];
chartlat="";
chartlat=BinaryRead[infile,"Real64"];
chartlong="";
chartlong=BinaryRead[infile,"Real64"];
chartmax="";
chartmax=BinaryRead[infile,"Real64"];
chartmin="";
chartmin=BinaryRead[infile,"Real64"];
chartTZ="";
chartTZ=BinaryRead[infile,"Integer16"];
chartsource="";
chartsource=ReadMyString[infile,10];
chartauthor="";
chartauthor=ReadMyString[infile,20];
chartlocalname="";
chartlocalname=ReadMyString[infile,20];
chartlocation="";
chartlocation=ReadMyString[infile,40];
chartchannels="";
chartchannels=BinaryRead[infile,"Integer16"];
chartnotelength="";
chartnotelength=BinaryRead[infile,"Integer32"];
chartnote="";
chartnote=ReadMyString[infile,chartnotelength];

(* Calculate parameters stored in header info and format data   *)
(* for display.                                                 *)
chartduration=Round2[(chartend-chartstart)*86400000]*0.001;
chartstartdate=DateList[(chartstart-2)*86400];
chartenddate=DateList[(chartend-2)*86400];
chartstartF=DateString[chartstartdate, "Day"]<>"-"<>DateString[chartstartdate, "MonthName"]<>"-"<>DateString[chartstartdate, "Year"]<>"  "<>DateString[chartstartdate, "Hour24"]<>":"<>DateString[chartstartdate, "Minute"]<>":"<>DateString[chartstartdate, "Second"]<>"."<>DateString[chartstartdate, "Millisecond"]<>" UTC";
chartendF=DateString[chartenddate, "Day"]<>"-"<>DateString[chartenddate, "MonthName"]<>"-"<>DateString[chartenddate, "Year"]<>"  "<>DateString[chartenddate, "Hour24"]<>":"<>DateString[chartenddate, "Minute"]<>":"<>DateString[chartenddate, "Second"]<>"."<>DateString[chartenddate, "Millisecond"]<>" UTC";
chartversionF=ToString[FromDigits[StringTake[chartversion,5]]]<>"."<>ToString[FromDigits[StringTake[chartversion,{6,7}]]]<>"."<>ToString[FromDigits[StringTake[chartversion,{8,10}]]];
freqlo=ToExpression[StringTake[chartnote,{ StringPosition[chartnote,"LOWF"][[1,2]]+1 ,StringPosition[chartnote,"LOWF"][[1,2]]+8}]]*10^-6;
freqhi=ToExpression[StringTake[chartnote,{ StringPosition[chartnote,"HIF"][[1,2]]+1 ,StringPosition[chartnote,"HIF"][[1,2]]+8}]]*10^-6;
channelsep=N[(freqhi-freqlo)/(chartchannels+1)*10^3]; (* channels+1 to account for RSS channel allocation bug *)
freqhi=N[freqhi-channelsep*10^-3]; (* account for RSS channel allocation bug *)
channeldwelltime=Round2[chartduration/(sweeps*chartchannels)*10^6];
sweeps=ToExpression[StringTake[chartnote,{ StringPosition[chartnote,"SWEEPS"][[1,2]]+1 ,StringPosition[chartnote,"ÿLOWF"][[1,1]]-1}]];
chartsweeprate=sweeps/chartduration;
If[StringPosition[chartsource,FromCharacterCode[0]][[1,1]]==1,chartsource="< not specified >"];
chartsweeps=Round2[chartsweeprate*plotduration];
shrinkratio=chartsweeps/chartwidth;
startoffset=Floor[86400*(chartstart-IntegerPart[chartstart])*chartsweeprate];

(* Read SPS data file into dataAraw and dataBraw matrices       *)
(* Iterate read function over sweeps and over the number of     *)
(* chart channels * 2, plus one EOR marker per sweep.           *)
(* NOTE: SPS data file ADC values are big endian (the byte      *)
(* order is reversed; i.e., MSB is the low byte)                *)
dataAraw=Table[0,{i,1,chartsweeps},{j,1,chartchannels}];
dataBraw=Table[0,{i,1,chartsweeps},{j,1,chartchannels}];
Do[
{
(* Read one input data record                                 *)
Do[
{
(* Read Side B data (side B comes first in the SPS data)    *)
dataBraw[[i,j]]=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
(* Read Side A data                                         *)
dataAraw[[i,j]]=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
}
(* cycle over 2 x number of channels *)
,{j,1,chartchannels}
];
(* Read EOR marker                                            *)
dummy=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
}
(* cycle over number of input sweeps to be averaged            *)
,{i,1+startoffset,sweeps+startoffset} 
];
Close[infile];


(*  RSS records channel 300 rolled over above channel 1.              *)
(*  These two lines shift channel 300 data back to where it belongs.  *)
dataAraw=RotateLeft[dataAraw,{0,1}];
dataBraw=RotateLeft[dataBraw,{0,1}];


(* Apply correction arrays to raw data                           *)
dataAcorr=dataAraw;
dataBcorr=dataBraw;
Do[dataAcorr[[i]]=dataAraw[[i]]+corrarray[[1]],{i,1,chartsweeps}];
Do[dataBcorr[[i]]=dataBraw[[i]]+corrarray[[2]],{i,1,chartsweeps}];


(* Do Stokes parameter subtractions on crrected data                    *)
dataAstokes=dataAcorr;
dataBstokes=dataBcorr;
Do[
{
A=dataAcorr[[i,j]];
B=dataBcorr[[i,j]];
Ax=0;
Bx=0;
(* Subtract LCP value from RCP value and store in RCP position *)
(* Multiply difference by 5 to make it easier to see in the spectrogram display *)
(* Keep values in 12-bit range of 0 to 4095 inclusive          *)
If[A-B≥0,If[5*(A-B)>4095,Ax=4095,Ax=5*(A-B)],Ax=0];
(* Subtract RCP value from LCP value and store in LCP position *)
If[B-A≥0,If[5*(B-A)>4095,Bx=4095,Bx=5*(B-A)],Bx=0];
(* Write Side A data                                      *)
dataAstokes[[i,j]]=Ax;
(* Write Side B data                                      *)
dataBstokes[[i,j]]=Bx;
}
,{i,1,chartsweeps},{j,1,chartchannels}];


(* Decimate raw data to fit within the spectrogram image width   *)
dataArawdecimated=Table[0,{i,1,chartwidth},{j,1,chartchannels}];
dataBrawdecimated=Table[0,{i,1,chartwidth},{j,1,chartchannels}];
Do[{
dataArawdecimated[[i,j]]=Mean[dataAraw[[Round2[1+(i-1)*shrinkratio];;Round2[i*shrinkratio],j]]];
dataBrawdecimated[[i,j]]=Mean[dataBraw[[Round2[1+(i-1)*shrinkratio];;Round2[i*shrinkratio],j]]];
},{j,1, chartchannels},{i,1, chartwidth}];


(* Decimate corrected data to fit within the spectrogram image width   *)
dataAcorrdecimated=Table[0,{i,1,chartwidth},{j,1,chartchannels}];
dataBcorrdecimated=Table[0,{i,1,chartwidth},{j,1,chartchannels}];
Do[{
dataAcorrdecimated[[i,j]]=Mean[dataAcorr[[Round2[1+(i-1)*shrinkratio];;Round2[i*shrinkratio],j]]];
dataBcorrdecimated[[i,j]]=Mean[dataBcorr[[Round2[1+(i-1)*shrinkratio];;Round2[i*shrinkratio],j]]];
},{j,1, chartchannels},{i,1, chartwidth}];


(* Decimate stokes data to fit within the spectrogram image width   *)
dataAstokesdecimated=Table[0,{i,1,chartwidth},{j,1,chartchannels}];
dataBstokesdecimated=Table[0,{i,1,chartwidth},{j,1,chartchannels}];
Do[{
dataAstokesdecimated[[i,j]]=Mean[dataAstokes[[Round2[1+(i-1)*shrinkratio];;Round2[i*shrinkratio],j]]];
dataBstokesdecimated[[i,j]]=Mean[dataBstokes[[Round2[1+(i-1)*shrinkratio];;Round2[i*shrinkratio],j]]];
},{j,1, chartchannels},{i,1, chartwidth}];


(* Apply color offset and color gain to raw data, just like RSS *) 
dataArawdecimatedadj=Round2[(dataArawdecimated-rawcoloroffset)*rawcolorgain];
dataBrawdecimatedadj=Round2[(dataBrawdecimated-rawcoloroffset)*rawcolorgain];
(* Coerce color-adjusted data to fit wihtin the n-bit range     *)
datalo=1;
datahi=2^nbits-1;  
(* Coerce Side A data                                          *)
Do[{
If[
dataArawdecimatedadj[[i,j]]≤datahi,
{If[
dataArawdecimatedadj[[i,j]]≥datalo,
{dataArawdecimatedadj[[i,j]]=dataArawdecimatedadj[[i,j]]},
{dataArawdecimatedadj[[i,j]]=datalo}
]},{dataArawdecimatedadj[[i,j]]=datahi}];
},{i,1,chartwidth},{j,1,chartchannels}];
(* Coerce Side B data                                          *)
Do[{
If[
dataBrawdecimatedadj[[i,j]]≤datahi,
{If[
dataBrawdecimatedadj[[i,j]]≥datalo,
{dataBrawdecimatedadj[[i,j]]=dataBrawdecimatedadj[[i,j]]},
{dataBrawdecimatedadj[[i,j]]=datalo}
]},{dataBrawdecimatedadj[[i,j]]=datahi}]
},{i,1,chartwidth},{j,1,chartchannels}];


(* Apply color offset and color gain to corrected data, just like RSS *) 
dataAcorrdecimatedadj=Round2[(dataAcorrdecimated-corrcoloroffset)*corrcolorgain];
dataBcorrdecimatedadj=Round2[(dataBcorrdecimated-corrcoloroffset)*corrcolorgain];
(* Coerce color-adjusted data to fit wihtin the n-bit range     *)
datalo=1;
datahi=2^nbits-1;  
(* Coerce Side A data                                          *)
Do[{
If[
dataAcorrdecimatedadj[[i,j]]≤datahi,
{If[
dataAcorrdecimatedadj[[i,j]]≥datalo,
{dataAcorrdecimatedadj[[i,j]]=dataAcorrdecimatedadj[[i,j]]},
{dataAcorrdecimatedadj[[i,j]]=datalo}
]},{dataAcorrdecimatedadj[[i,j]]=datahi}];
},{i,1,chartwidth},{j,1,chartchannels}];
(* Coerce Side B data                                          *)
Do[{
If[
dataBcorrdecimatedadj[[i,j]]≤datahi,
{If[
dataBcorrdecimatedadj[[i,j]]≥datalo,
{dataBcorrdecimatedadj[[i,j]]=dataBcorrdecimatedadj[[i,j]]},
{dataBcorrdecimatedadj[[i,j]]=datalo}
]},{dataBcorrdecimatedadj[[i,j]]=datahi}]
},{i,1,chartwidth},{j,1,chartchannels}];


(* Apply color offset and color gain to stokes data, just like RSS *) 
dataAstokesdecimatedadj=Round2[(dataAstokesdecimated-stokescoloroffset)*stokescolorgain];
dataBstokesdecimatedadj=Round2[(dataBstokesdecimated-stokescoloroffset)*stokescolorgain];
(* Coerce color-adjusted data to fit wihtin the n-bit range     *)
datalo=1;
datahi=2^nbits-1;  
(* Coerce Side A data                                          *)
Do[{
If[
dataAstokesdecimatedadj[[i,j]]≤datahi,
{If[
dataAstokesdecimatedadj[[i,j]]≥datalo,
{dataAstokesdecimatedadj[[i,j]]=dataAstokesdecimatedadj[[i,j]]},
{dataAstokesdecimatedadj[[i,j]]=datalo}
]},{dataAstokesdecimatedadj[[i,j]]=datahi}];
},{i,1,chartwidth},{j,1,chartchannels}];
(* Coerce Side B data                                          *)
Do[{
If[
dataBstokesdecimatedadj[[i,j]]≤datahi,
{If[
dataBstokesdecimatedadj[[i,j]]≥datalo,
{dataBstokesdecimatedadj[[i,j]]=dataBstokesdecimatedadj[[i,j]]},
{dataBstokesdecimatedadj[[i,j]]=datalo}
]},{dataBstokesdecimatedadj[[i,j]]=datahi}]
},{i,1,chartwidth},{j,1,chartchannels}];


plot1=ArrayPlot[Transpose[dataArawdecimatedadj],
PixelConstrained->True,
ImageSize->{chartwidth,chartchannels},
ColorFunctionScaling->False,
ColorFunction->Function[{z},RGBColor[ColorsInT[[z]]]],
Background->None,
Frame->False
];

plot2=ArrayPlot[Transpose[dataBrawdecimatedadj],
PixelConstrained->True,
ImageSize->{chartwidth,chartchannels},
ColorFunctionScaling->False,
ColorFunction->Function[{z},RGBColor[ColorsInT[[z]]]],
Background->None,
Frame->False
];

plot3=ArrayPlot[Transpose[dataAcorrdecimatedadj],
PixelConstrained->True,
ImageSize->{chartwidth,chartchannels},
ColorFunctionScaling->False,
ColorFunction->Function[{z},RGBColor[ColorsInT[[z]]]],
Background->None,
Frame->False
];

plot4=ArrayPlot[Transpose[dataBcorrdecimatedadj],
PixelConstrained->True,
ImageSize->{chartwidth,chartchannels},
ColorFunctionScaling->False,
ColorFunction->Function[{z},RGBColor[ColorsInT[[z]]]],
Background->None,
Frame->False
];

plot5=ArrayPlot[Transpose[dataAstokesdecimatedadj],
PixelConstrained->True,
ImageSize->{chartwidth,chartchannels},
ColorFunctionScaling->False,
ColorFunction->Function[{z},RGBColor[ColorsInT[[z]]]],
Background->None,
Frame->False
];

plot6=ArrayPlot[Transpose[dataBstokesdecimatedadj],
PixelConstrained->True,
ImageSize->{chartwidth,chartchannels},
ColorFunctionScaling->False,
ColorFunction->Function[{z},RGBColor[ColorsInT[[z]]]],
Background->None,
Frame->False
];

(* export PNG files *)
Export[spgmAfilename,plot1];
Export[spgmBfilename,plot2];
Export[spgmAfilenamecorr,plot3];
Export[spgmBfilenamecorr,plot4];
Export[spgmAfilenamecorrstokes,plot5];
Export[spgmBfilenamecorrstokes,plot6];


},{z,1,filecount}
];

(* Print ens time stamp *)
Print["Processing Complete on "<>DateString[]]

chartstart

chartstart*86400

chartend*86400

(chartend-chartstart)*86400

86400*(chartstart-IntegerPart[chartstart])

86400*(chartend-IntegerPart[chartend])

Files

Floor[86400*(chartstart-IntegerPart[chartstart])*chartsweeprate]

Floor