Dual Polarization Spectrograph Data File Utilities
Dave Typinski, AJ4CO Observatory
May, 2013

DPS Side Swapper - Swaps Side A data with Side B data and writes new SPS file

(********************************************************************************************)
(*                                                                                          *)
(*  Swap Side A with Side B data and write new SPS file                                     *)
(*                                                                                          *)
(*                                                                                          *)
(********************************************************************************************)

(* Open one files for input                                     *)
filename="R:\\Observation Records\\AJ4CO_DPS-XTFD\\1305\\130503\\130503171500.sps";
infile=OpenRead[filename,BinaryFormat->True];
outfile=OpenWrite[StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>"_rev_pol.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
];

(* 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]

DPS Channel Rotator - Rotates channel 200 back down below channel 199 and writes new SPS file

(********************************************************************************************)
(*                                                                                          *)
(*  RSS and DPS data have channel 200 rolled over above channel 1.                          *)
(*  This routine puts channel 200 data back where it belongs.                               *)
(*                                                                                          *)
(********************************************************************************************)

(* Open one file for input                                      *)
filename="R:\\Observation Records\\AJ4CO_DPS-XTFD\\1305\\130503\\130503085000.sps";
infile=OpenRead[filename,BinaryFormat->True];
outfile=OpenWrite[StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>"_r200.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
];

(* Read 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}]];

(* Create sweep arrays, then populate, rotate, 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 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];
Close[outfile];
Print["Done."]

DPS File Concatenator - Concatenates two DPS SPS files

(* File names                                                   *)
filename="R:\\Observation Records\\AJ4CO Observatory\\DPS\\1409\\140912\\140912000000.sps";
filename2="R:\\Observation Records\\AJ4CO Observatory\\DPS\\1409\\140912\\140912195222.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 Continuum Inspector - Plots Side A and Side B average continua of entire input file, channel matrix list, saves report as PDF

(********************************************************************************************)
(*                                                                                          *)
(*  Shows continuum of channel averages over SHORT input file duration                      *)
(*  This routine doe NOT perform any decimation                                             *)
(*                                                                                          *)
(********************************************************************************************)

(* Open one file for input                                      *)
filename="R:\\Observation Records\\AJ4CO_DPS-XTFD\\1305\\130531\\130531222140.sps";
infile=OpenRead[filename,BinaryFormat->True];
outfilename=StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>"_continua_inspection.pdf";
outfilenamecsv=StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>"_continua_inspection.csv";

(* 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
];

(* 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;

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

(* Extract sweep range from note entry & calc channel separation *)
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=(freqhi-freqlo)/chartchannels*10^3;

(* Create sweep arrays, then populate, rotate, write to file    *)
(* iterating over all sweeps in source file.                    *)
indataB=Table[0,{i,1,chartchannels}];
indataA=Table[0,{i,1,chartchannels}];
avgB=Table[0,{i,1,chartchannels}];
avgA=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}
];
(* Add sweep to average                                      *)
avgB=avgB+indataB;
avgA=avgA+indataA;
(* Read EOR marker                                           *)
eor=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
}
,{j,1,sweeps}
];
Close[infile];

(* Average the data                                             *)
avgB=Round2[avgB/sweeps];
avgA=Round2[avgA/sweeps];

(*  RSS and DPS data have channel 200 rolled over above channel 1. *)
(*  These two lines put channel 200 data back where it belongs.    *)
avgB=RotateLeft[avgB];
avgA=RotateLeft[avgA];

(* Create channel frequency list                                *)
dataFM=Table[ToString[PaddedForm[N[freqhi-i*channelsep*10^-3],{5,3},NumberPadding->{"","0"}]],{i,1,chartchannels}];

(* Create a continuum plot                                      *)
charttitle=Column[{"",
Style["RCP & LCP Continua – average of "<>ToString[sweeps]<>" sweeps ("<>ToString[Round2[chartduration]]<>" seconds of observation data)",Plain,FontSize->12,FontFamily->Helvetica],
Row[{
Style["RCP = Green",Plain,FontSize->12,FontFamily->Helvetica,RGBColor[0,150/255,100/255]],
Style["LCP = RED",Plain,FontSize->12,FontFamily->Helvetica,RGBColor[190/255,39/255,39/255]]
},"        "]
},Center];
plot1=ListLinePlot[{avgA,avgB},
DataRange->{freqlo,freqhi},
ScalingFunctions->{"Reverse"},
PlotRange->{1500,4000},
PlotStyle->{RGBColor[0,150/255,100/255],RGBColor[190/255,39/255,39/255]},
PlotMarkers->None,
Axes->None,
GridLines->None,
Frame->True,
RotateLabel-> True,
FrameStyle->{{Directive[Thin,FontFamily->Helvetica,11],Thin},{Directive[Thin,FontFamily->Helvetica,11],Thin}},FrameLabel->{Style["Frequency (MHz)",FontSize->11,FontFamily->Helvetica],Style["Average ADC Value (mV)",FontSize->11,FontFamily->Helvetica],charttitle,""}
]

(* Prepend column labeles and position numbers                  *)
A=Prepend[avgA,"  Side A  "];
B=Prepend[avgB,"  Side B  "];
Frqs=Prepend[dataFM,"  Freq (MHz)  "];
Ch=Prepend[Table[i,{i,1,chartchannels}],"Channel"];
Dummy=Table["        ",{i,1,chartchannels+1}];

(* Create output array                                          *)
output=Table["",{i,1,50}];
output[[1]]=Style["DPS Average Continuum Inspector",Larger,Darker[Darker[Blue]],Bold];
output[[2]]="Dave Typinski, May 2013";
output[[10]]=Style["Analysis",Bold];
output[[12]]="         Input file: "<> filename;
output[[13]]="      Data duration: "<>ToString[PaddedForm[chartduration,{10,3},NumberPadding->{"","0"}]]<>" s";
output[[14]]="             Sweeps: "<>ToString[sweeps];
output[[16]]=Style["  NOTE: N-S feed line on 6ft-6in (90° @ 25 MHz) phasing line on CP test jig.",Bold];
output[[20]]=Rasterize[plot1,ImageSize->650];
output[[22]]=TableForm[Transpose[{Dummy,Ch,Frqs,A,B}],TableAlignments->Center];

(* Output the result                                            *)
g=Panel[Column[{Column[output[[10;;11]]],Column[output[[12;;17]],Alignment->":"],output[[20]],output[[21]],output[[22]]}],Column[output[[1;;3]]]];
k=Column[{g}];
g
Export[outfilename,g]

outcsvdata=Transpose[{Ch,Frqs,A,B}];
outcsvdata=Prepend[outcsvdata,"Values below averaged from "<>ToString[sweeps]<>" sweeps ("<>ToString[PaddedForm[chartduration,{10,3},NumberPadding->{"","0"}]]<>" seconds of observation)"];
outcsvdata=Prepend[outcsvdata,"Frequency sweep range "<>ToString[freqlo]<>" to "<>ToString[freqhi]<>" MHz"];
outcsvdata=Prepend[outcsvdata,"Continuum data from "<>filename];
outcsvdata=Prepend[outcsvdata,"CSV file created "<>DateString[{"Day","-","MonthNameShort","-","Year"," ","Hour24",":","Minute"," Local"}]];
outcsvdata=Prepend[outcsvdata,outfilenamecsv];
outcsvdata=Prepend[outcsvdata,"DPS Continuum Inspector Report"];
Export[outfilenamecsv,outcsvdata]

DPS Step Calibration Inspector - Plots Side A and Side B data from multiple input files, saves report as PDF & data as CSV

(********************************************************************************************)
(*                                                                                          *)
(*  Shows continuum of channel averages over input file duration for multiple files         *)
(*                                                                                          *)
(*                                                                                          *)
(********************************************************************************************)

(* Set directory for input                                      *)
filebase="R:\\Observation Records\\AJ4CO_DPS-XTFD\\1305\\130525\\";

(* Set number of SPS files to analyze                           *)
(* Files must be named 00.sps, 01.sps, 02.sps, etc.             *)
filecount=21;

(* Create output filenames                                      *)
rptfilename=filebase<>"DPS Step cal report ";
csvfilename=filebase<>"DPS Step cal data.csv";

(* Create empty data matrices                                   *)
avgAmatrix=Table[0,{i,1,filecount}];
avgBmatrix=Table[0,{i,1,filecount}];

(* Create empty panel array                                     *)
panels=Table[0,{i,1,filecount}];

(* Create empty plots arrays                                    *)
plots=Table[0,{i,1,filecount}];
plotsR=Table[0,{i,1,filecount}];
plotsL=Table[0,{i,1,filecount}];

(* 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]];

(* Main loop to cycle through all input files                   *)
Do[{

(* Create filename and open for read                            *)
filename=filebase<>IntegerString[k-1,10,2]<>".sps";
infile=OpenRead[filename,BinaryFormat->True];

(* 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];
chartstartF=DateString[chartstartdate, "Day"]<>"-"<>DateString[chartstartdate, "MonthName"]<>"-"<>DateString[chartstartdate, "Year"];

(* Read number of frequency sweeps from neader notes            *)
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=(freqhi-freqlo)/chartchannels*10^3;
sweeps=ToExpression[StringTake[chartnote,{ StringPosition[chartnote,"SWEEPS"][[1,2]]+1 ,StringPosition[chartnote,"ÿLOWF"][[1,1]]-1}]];

(* Create sweep arrays, then populate, rotate, write to file    *)
(* iterating over all sweeps in source file.                    *)
indataB=Table[0,{i,1,chartchannels}];
indataA=Table[0,{i,1,chartchannels}];
avgB=Table[0,{i,1,chartchannels}];
avgA=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}
];
(* Add sweep to average                                      *)
avgB=avgB+indataB;
avgA=avgA+indataA;
(* Read EOR marker                                           *)
eor=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
}
,{j,1,sweeps}
];
Close[infile];

(* Average the data                                             *)
avgB=Round2[avgB/sweeps];
avgA=Round2[avgA/sweeps];

(*  RSS and DPS data have channel 200 rolled over above channel 1. *)
(*  These two lines put channel 200 data back where it belongs.    *)

avgB=RotateLeft[avgB];
avgA=RotateLeft[avgA];

avgAmatrix[[k]]=avgA;
avgBmatrix[[k]]=avgB;

(* Create a continuum plot                                      *)
charttitle="";
plots[[k]]=ListLinePlot[{avgA,avgB},
DataRange->{1,chartchannels},
PlotRange->{1000,4100},
PlotStyle->{RGBColor[0,150/255,100/255],RGBColor[190/255,39/255,39/255]},
PlotMarkers->None,
Axes->None,
GridLines->None,
Frame->True,
RotateLabel-> True,
FrameStyle->{{Directive[Thin,FontFamily->Helvetica,11],Thin},{Directive[Thin,FontFamily->Helvetica,11],Thin}},FrameLabel->{Style["Channel Number",FontSize->11,FontFamily->Helvetica],Style["Average ADC Value (mV)",FontSize->11,FontFamily->Helvetica],charttitle,""}
];

plotsR[[k]]=ListLinePlot[{avgA},
DataRange->{1,chartchannels},
PlotRange->{1000,4100},
PlotStyle->{RGBColor[0,150/255,100/255]},
PlotMarkers->None,
Axes->None,
GridLines->None,
Frame->True,
RotateLabel-> True,
FrameStyle->{{Directive[Thin,FontFamily->Helvetica,11],Thin},{Directive[Thin,FontFamily->Helvetica,11],Thin}},FrameLabel->{Style["Channel Number",FontSize->11,FontFamily->Helvetica],Style["Average ADC Value (mV)",FontSize->11,FontFamily->Helvetica],charttitle,""}
];

plotsL[[k]]=ListLinePlot[{avgB},
DataRange->{1,chartchannels},
PlotRange->{1000,4100},
PlotStyle->{RGBColor[190/255,39/255,39/255]},
PlotMarkers->None,
Axes->None,
GridLines->None,
Frame->True,
RotateLabel-> True,
FrameStyle->{{Directive[Thin,FontFamily->Helvetica,11],Thin},{Directive[Thin,FontFamily->Helvetica,11],Thin}},FrameLabel->{Style["Channel Number",FontSize->11,FontFamily->Helvetica],Style["Average ADC Value (mV)",FontSize->11,FontFamily->Helvetica],charttitle,""}
];


(* Create output array                                          *)
output=Table["",{i,1,50}];
output[[1]]=Style["DPS Step Calibration Inspector",Larger,Darker[Darker[Blue]],Bold];
output[[2]]="Dave Typinski, May 2013";
output[[10]]=Style["Analysis",Bold];
output[[12]]="         Input file: "<> filename;
output[[13]]="      Data duration: "<>ToString[PaddedForm[chartduration,{10,3},NumberPadding->{"","0"}]]<>" s";
output[[14]]="             Sweeps: "<>ToString[sweeps];
output[[16]]="Continuum Plot: Average of "<>ToString[sweeps]<>" Sweeps, "<>
Which[
k≤18,IntegerString[3*(k-1)]<>" dB Attenuation",
k==19,"00 dB Attenuation (drift check)",
k==20,"50 ohm term on splitter input",
k==21,"50 ohm term on each DPS input"
];
output[[18]]=Row[{
Style["              Green = RCP",RGBColor[0,150/255,100/255]],
"       ",
Style["Red = LCP",RGBColor[190/255,39/255,39/255]]
}];

output[[20]]=Rasterize[plots[[k]],ImageSize->650];

(* Output the result                                            *)
panels[[k]]=Panel[Column[{
Column[output[[10;;11]]],
Column[output[[12;;17]],Alignment->":"],
Column[output[[18;;19]]],
output[[20]],
"",
Style["RCP channels, averaged data values",RGBColor[0,150/255,100/255]],
avgA,
"",
Style["LCP channels, averaged data values",RGBColor[190/255,39/255,39/255]],
avgB
}],
Column[output[[1;;3]]]];

Export[rptfilename<>IntegerString[k-1,10,2]<>".pdf",panels[[k]]];

},{k,1,filecount}
];

(* Output the result                                            *)

g=Show[Table[plots[[i]],{i,1,17,2}]];
output[[20]]=Rasterize[g,ImageSize->650];
panelmix=Panel[Column[{
Column[output[[10;;11]]],
"Combined continua plot, 0 dB to 48 dB attenuation, 6 dB steps",
Column[output[[17;;19]]],
output[[20]]
}],
Column[output[[1;;3]]]];
Export[rptfilename<>"combined plot"<>".pdf",panelmix];

g=Show[Table[plotsR[[i]],{i,1,17,2}]];
output[[20]]=Rasterize[g,ImageSize->650];
panelmix=Panel[Column[{
Column[output[[10;;11]]],
"Combined continua plot, 0 dB to 48 dB attenuation, 6 dB steps, RCP only",
Column[output[[17;;19]]],
output[[20]]
}],
Column[output[[1;;3]]]];
Export[rptfilename<>"combined plot RCP"<>".pdf",panelmix];

g=Show[Table[plotsL[[i]],{i,1,17,2}]];
output[[20]]=Rasterize[g,ImageSize->650];
panelmix=Panel[Column[{
Column[output[[10;;11]]],
"Combined continua plot, 0 dB to 48 dB attenuation, 6 dB steps, LCP only",
Column[output[[17;;19]]],
output[[20]]
}],
Column[output[[1;;3]]]];
Export[rptfilename<>"combined plot LCP"<>".pdf",panelmix];

lplots=Table[0,{i,1,5}];

Do[
lplots[[m]]=
ListLinePlot[{avgAmatrix[[1;;18,Which[m-1<1,1,m-1>0,50*(m-1)]]]},
DataRange->{0,51},
PlotRange->{1000,4100},
PlotStyle->{RGBColor[0,150/255,100/255]},
PlotMarkers->None,
Axes->None,
GridLines->None,
Frame->True,
RotateLabel-> True,
FrameStyle->{{Directive[Thin,FontFamily->Helvetica,11],Thin},{Directive[Thin,FontFamily->Helvetica,11],Thin}},FrameLabel->{Style["Attenuation (dB)",FontSize->11,FontFamily->Helvetica],Style["Average ADC Value (mV)",FontSize->11,FontFamily->Helvetica],charttitle,""}
],
{m,1,5}
];
g=Show[Table[lplots[[i]],{i,1,5}]];
output[[20]]=Rasterize[g,ImageSize->650];
panelmix=Panel[Column[{
Column[output[[10;;11]]],
"Response plot (dynamic range), channels 1, 50, 100, 150, and 200 – RCP only",
Column[output[[17;;19]]],
output[[20]],
"",
Style["RCP channels 1, 50, 100, 150, and 200, averaged data values",RGBColor[0,150/255,100/255]],
MatrixForm[
Prepend[
MapThread[Prepend,{Table[avgAmatrix[[1;;18,Which[v-1<1,1,v-1>0,50*(v-1)]]],{v,1,5}],{"Channel 001:","Channel 050:","Channel 100:","Channel 150:","Channel 200:"}}],
{"            "," 0dB"," 3dB"," 6dB"," 9dB","12dB","15dB","18dB","21dB","24dB","27dB","30dB","33dB","36dB","39dB","42dB","45dB","48dB","51dB"}
]
]
}],
Column[output[[1;;3]]]];
Export[rptfilename<>"dyn range plot RCP"<>".pdf",panelmix];

Do[
lplots[[m]]=
ListLinePlot[{avgBmatrix[[1;;18,Which[m-1<1,1,m-1>0,50*(m-1)]]]},
DataRange->{0,51},
PlotRange->{1000,4100},
PlotStyle->{RGBColor[190/255,39/255,39/255]},
PlotMarkers->None,
Axes->None,
GridLines->None,
Frame->True,
RotateLabel-> True,
FrameStyle->{{Directive[Thin,FontFamily->Helvetica,11],Thin},{Directive[Thin,FontFamily->Helvetica,11],Thin}},FrameLabel->{Style["Attenuation (dB)",FontSize->11,FontFamily->Helvetica],Style["Average ADC Value (mV)",FontSize->11,FontFamily->Helvetica],charttitle,""}
],
{m,1,5}
];
g=Show[Table[lplots[[i]],{i,1,5}]];
output[[20]]=Rasterize[g,ImageSize->650];
panelmix=Panel[Column[{
Column[output[[10;;11]]],
"Response plot (dynamic range), channels 1, 50, 100, 150, and 200 – LCP only",
Column[output[[17;;19]]],
output[[20]],
"",
Style["RCP channels 1, 50, 100, 150, and 200, averaged data values",RGBColor[0,150/255,100/255]],
MatrixForm[
Prepend[
MapThread[Prepend,{Table[avgBmatrix[[1;;18,Which[v-1<1,1,v-1>0,50*(v-1)]]],{v,1,5}],{"Channel 001:","Channel 050:","Channel 100:","Channel 150:","Channel 200:"}}],
{"            "," 0dB"," 3dB"," 6dB"," 9dB","12dB","15dB","18dB","21dB","24dB","27dB","30dB","33dB","36dB","39dB","42dB","45dB","48dB","51dB"}
]
]
}],
Column[output[[1;;3]]]];
Export[rptfilename<>"dyn range plot LCP"<>".pdf",panelmix];

Do[
lplots[[m]]=
ListLinePlot[{avgAmatrix[[1;;18,Which[m-1<1,1,m-1>0,50*(m-1)]]],avgBmatrix[[1;;18,Which[m-1==0,1,m-1>0,50*(m-1)]]]},
DataRange->{0,51},
PlotRange->{1000,4100},
PlotStyle->{RGBColor[0,150/255,100/255],RGBColor[190/255,39/255,39/255]},
PlotMarkers->None,
Axes->None,
GridLines->None,
Frame->True,
RotateLabel-> True,
FrameStyle->{{Directive[Thin,FontFamily->Helvetica,11],Thin},{Directive[Thin,FontFamily->Helvetica,11],Thin}},FrameLabel->{Style["Attenuation (dB)",FontSize->11,FontFamily->Helvetica],Style["Average ADC Value (mV)",FontSize->11,FontFamily->Helvetica],charttitle,""}
],
{m,1,5}
];
g=Show[Table[lplots[[i]],{i,1,5}]];
output[[20]]=Rasterize[g,ImageSize->650];
panelmix=Panel[Column[{
Column[output[[10;;11]]],
"Response plot (dynamic range), channels 1, 50, 100, 150, and 200 – RCP and LCP",
Column[output[[17;;19]]],
output[[20]],
"",
Style["RCP channels 1, 50, 100, 150, and 200, averaged data values",RGBColor[0,150/255,100/255]],
MatrixForm[
Prepend[
MapThread[Prepend,{Table[avgAmatrix[[1;;18,Which[v-1<1,1,v-1>0,50*(v-1)]]],{v,1,5}],{"Channel 001:","Channel 050:","Channel 100:","Channel 150:","Channel 200:"}}],
{"            "," 0dB"," 3dB"," 6dB"," 9dB","12dB","15dB","18dB","21dB","24dB","27dB","30dB","33dB","36dB","39dB","42dB","45dB","48dB","51dB"}
]
],
"",
Style["LCP channels 1, 50, 100, 150, and 200, averaged data values",RGBColor[190/255,39/255,39/255]],
MatrixForm[
Prepend[
MapThread[Prepend,{Table[avgBmatrix[[1;;18,Which[v-1<1,1,v-1>0,50*(v-1)]]],{v,1,5}],{"Channel 001:","Channel 050:","Channel 100:","Channel 150:","Channel 200:"}}],
{"            "," 0dB"," 3dB"," 6dB"," 9dB","12dB","15dB","18dB","21dB","24dB","27dB","30dB","33dB","36dB","39dB","42dB","45dB","48dB","51dB"}
]
]
}],
Column[output[[1;;3]]]];
Export[rptfilename<>"dyn range plot"<>".pdf",panelmix];

(* Average mV per dB analysis routines *)
respA=Table[0,{i,1,5}];
respB=Table[0,{i,1,5}];
Do[{
respA[[m]]=
avgAmatrix[[1;;18,Which[m-1<1,1,m-1>0,50*(m-1)]]];
respB[[m]]=
avgBmatrix[[1;;18,Which[m-1<1,1,m-1>0,50*(m-1)]]]
},
{m,1,5}
];
respAc=Table[0,{i,1,5},{j,1,17}];
respBc=Table[0,{i,1,5},{j,1,17}];
Do[{
Do[{
respAc[[v,m]]=Round2[(respA[[v,m]]-respA[[v,m+1]])/3];
respBc[[v,m]]=Round2[(respB[[v,m]]-respB[[v,m+1]])/3];
},
{m,1,17}
];
},
{v,1,5}
];

Do[{
lplots[[m]]=
ListLinePlot[{respAc[[m]]},
DataRange->{1.5,49.5},
PlotRange->{0,60},
PlotStyle->{RGBColor[0,150/255,100/255]},
PlotMarkers->None,
Axes->None,
GridLines->None,
Frame->True,
RotateLabel-> True,
FrameStyle->{{Directive[Thin,FontFamily->Helvetica,11],Thin},{Directive[Thin,FontFamily->Helvetica,11],Thin}},FrameLabel->{Style["Attenuation (dB)",FontSize->11,FontFamily->Helvetica],Style["ADC mV per dB",FontSize->11,FontFamily->Helvetica],charttitle,""}
];
},
{m,1,5}
];
g=Show[Table[lplots[[i]],{i,1,5}]];
output[[20]]=Rasterize[g,ImageSize->650];
panelmix=Panel[Column[{
Column[output[[10;;11]]],
"Detector response in ADC mV per dB, channels 1, 50, 100, 150, and 200 – RCP only",
Column[output[[17;;19]]],
output[[20]],
"",
Style["RCP channels 1, 50, 100, 150, and 200, detector response data values in mV/dB",RGBColor[0,150/255,100/255]],
MatrixForm[
Prepend[
MapThread[Prepend,{respAc,{"Ch 001:","Ch 050:","Ch 100:","Ch 150:","Ch 200:"}}],
{"            ","1.5dB","4.5dB","7.5dB","11.5dB","13.5dB","16.5dB","19.5dB","22.5dB","25.5dB","28.5dB","31.5dB","34.5dB","37.5dB","40.5dB","43.5dB","46.5dB","49.5dB"}
]
],
"",
Style["RCP channels 1, 50, 100, 150, and 200, averaged data values in mV",RGBColor[0,150/255,100/255]],
MatrixForm[
Prepend[
MapThread[Prepend,{Table[avgAmatrix[[1;;18,Which[v-1<1,1,v-1>0,50*(v-1)]]],{v,1,5}],{"Channel 001:","Channel 050:","Channel 100:","Channel 150:","Channel 200:"}}],
{"            "," 0dB"," 3dB"," 6dB"," 9dB","12dB","15dB","18dB","21dB","24dB","27dB","30dB","33dB","36dB","39dB","42dB","45dB","48dB","51dB"}
]
]
}],
Column[output[[1;;3]]]];
Export[rptfilename<>"detector response plot RCP"<>".pdf",panelmix];

Do[{
lplots[[m]]=
ListLinePlot[{respBc[[m]]},
DataRange->{1.5,49.5},
PlotRange->{0,60},
PlotStyle->{RGBColor[190/255,39/255,39/255]},
PlotMarkers->None,
Axes->None,
GridLines->None,
Frame->True,
RotateLabel-> True,
FrameStyle->{{Directive[Thin,FontFamily->Helvetica,11],Thin},{Directive[Thin,FontFamily->Helvetica,11],Thin}},FrameLabel->{Style["Attenuation (dB)",FontSize->11,FontFamily->Helvetica],Style["ADC mV per dB",FontSize->11,FontFamily->Helvetica],charttitle,""}
];
},
{m,1,5}
];
g=Show[Table[lplots[[i]],{i,1,5}]];
output[[20]]=Rasterize[g,ImageSize->650];
panelmix=Panel[Column[{
Column[output[[10;;11]]],
"Detector response in ADC mV per dB, channels 1, 50, 100, 150, and 200 – LCP only",
Column[output[[17;;19]]],
output[[20]],
"",
Style["LCP channels 1, 50, 100, 150, and 200, detector response data values in mV/dB",RGBColor[190/255,39/255,39/255]],
MatrixForm[
Prepend[
MapThread[Prepend,{respBc,{"Ch 001:","Ch 050:","Ch 100:","Ch 150:","Ch 200:"}}],
{"            ","1.5dB","4.5dB","7.5dB","11.5dB","13.5dB","16.5dB","19.5dB","22.5dB","25.5dB","28.5dB","31.5dB","34.5dB","37.5dB","40.5dB","43.5dB","46.5dB","49.5dB"}
]
],
"",
Style["LCP channels 1, 50, 100, 150, and 200, averaged data values in mV",RGBColor[190/255,39/255,39/255]],
MatrixForm[
Prepend[
MapThread[Prepend,{Table[avgBmatrix[[1;;18,Which[v-1<1,1,v-1>0,50*(v-1)]]],{v,1,5}],{"Channel 001:","Channel 050:","Channel 100:","Channel 150:","Channel 200:"}}],
{"            "," 0dB"," 3dB"," 6dB"," 9dB","12dB","15dB","18dB","21dB","24dB","27dB","30dB","33dB","36dB","39dB","42dB","45dB","48dB","51dB"}
]
]
}],
Column[output[[1;;3]]]];
Export[rptfilename<>"detector response plot LCP"<>".pdf",panelmix];

avgrespA=Round2[Mean[respAc]];
avgrespB=Round2[Mean[respBc]];

w=ListLinePlot[{avgrespA,avgrespB},
DataRange->{1.5,49.5},
PlotRange->{0,60},
PlotStyle->{RGBColor[0,150/255,100/255],RGBColor[190/255,39/255,39/255]},
PlotMarkers->None,
Axes->None,
GridLines->None,
Frame->True,
RotateLabel-> True,
FrameStyle->{{Directive[Thin,FontFamily->Helvetica,11],Thin},{Directive[Thin,FontFamily->Helvetica,11],Thin}},FrameLabel->{Style["Attenuation (dB)",FontSize->11,FontFamily->Helvetica],Style["ADC mV per dB",FontSize->11,FontFamily->Helvetica],charttitle,""}
];
output[[20]]=Rasterize[w,ImageSize->650];
panelmix=Panel[Column[{
Column[output[[10;;11]]],
"Average detector response in ADC mV per dB, channels 1, 50, 100, 150, and 200",
Column[output[[17;;19]]],
output[[20]],
"",
"Channels 1, 50, 100, 150, and 200, average detector response data values in mV/dB",
MatrixForm[
Prepend[
MapThread[Prepend,{{avgrespA,avgrespB},{Style["Average RCP",RGBColor[0,150/255,100/255]],Style["Average LCP",RGBColor[190/255,39/255,39/255]]}}],
{"            ","1.5dB","4.5dB","7.5dB","11.5dB","13.5dB","16.5dB","19.5dB","22.5dB","25.5dB","28.5dB","31.5dB","34.5dB","37.5dB","40.5dB","43.5dB","46.5dB","49.5dB"}
]
],
"",
Style["RCP channels 1, 50, 100, 150, and 200, detector response data values in mV/dB",RGBColor[0,150/255,100/255]],
MatrixForm[
Prepend[
MapThread[Prepend,{respAc,{"Ch 001:","Ch 050:","Ch 100:","Ch 150:","Ch 200:"}}],
{"            ","1.5dB","4.5dB","7.5dB","11.5dB","13.5dB","16.5dB","19.5dB","22.5dB","25.5dB","28.5dB","31.5dB","34.5dB","37.5dB","40.5dB","43.5dB","46.5dB","49.5dB"}
]
],
"",
Style["RCP channels 1, 50, 100, 150, and 200, averaged data values in mV",RGBColor[0,150/255,100/255]],
MatrixForm[
Prepend[
MapThread[Prepend,{Table[avgAmatrix[[1;;18,Which[v-1<1,1,v-1>0,50*(v-1)]]],{v,1,5}],{"Channel 001:","Channel 050:","Channel 100:","Channel 150:","Channel 200:"}}],
{"            "," 0dB"," 3dB"," 6dB"," 9dB","12dB","15dB","18dB","21dB","24dB","27dB","30dB","33dB","36dB","39dB","42dB","45dB","48dB","51dB"}
]
],
"",
Style["LCP channels 1, 50, 100, 150, and 200, detector response data values in mV/dB",RGBColor[190/255,39/255,39/255]],
MatrixForm[
Prepend[
MapThread[Prepend,{respBc,{"Ch 001:","Ch 050:","Ch 100:","Ch 150:","Ch 200:"}}],
{"            ","1.5dB","4.5dB","7.5dB","11.5dB","13.5dB","16.5dB","19.5dB","22.5dB","25.5dB","28.5dB","31.5dB","34.5dB","37.5dB","40.5dB","43.5dB","46.5dB","49.5dB"}
]
],
"",
Style["LCP channels 1, 50, 100, 150, and 200, averaged data values in mV",RGBColor[190/255,39/255,39/255]],
MatrixForm[
Prepend[
MapThread[Prepend,{Table[avgBmatrix[[1;;18,Which[v-1<1,1,v-1>0,50*(v-1)]]],{v,1,5}],{"Channel 001:","Channel 050:","Channel 100:","Channel 150:","Channel 200:"}}],
{"            "," 0dB"," 3dB"," 6dB"," 9dB","12dB","15dB","18dB","21dB","24dB","27dB","30dB","33dB","36dB","39dB","42dB","45dB","48dB","51dB"}
]
]
}],
Column[output[[1;;3]]]];
Export[rptfilename<>"detector response plot"<>".pdf",panelmix];


(* Create channel frequency list                                *)
datafc=Table[ToString[PaddedForm[N[freqhi-i*channelsep*10^-3],{5,3},NumberPadding->{"","0"}]],{i,1,chartchannels}];
datafc=Prepend[datafc,"f_c (MHz)"];
datafc=Prepend[datafc,""];
datafc=Prepend[datafc,""];

(* Create channel number list *)
datacn=Table[i,{i,1,chartchannels}];
datacn=Prepend[datacn,"Channel"];
datacn=Prepend[datacn,"Temp @ DPS inputs (kK):"];
datacn=Prepend[datacn,"Attenuator setting (dB):"];

Amat=avgAmatrix;
AmatT=Transpose[Amat];
AmatT=Prepend[AmatT,Table["RCP mV",{i,1,21}]];
AmatT=Prepend[AmatT,{143800,72080,36120,18110,9075,4549,2280,1144,573.7,288.1,145.0,73.23,37.28,19.26,10.23,5.702,3.434,2.297,143800,0.290,0.290}];
AmatT=Prepend[AmatT,{0,3,6,9,12,15,18,21,24,27,30,33,36,39,42,45,48,51,0,"50 ohms on MC","Kay @ 101 dB"}];
AmatTT=Transpose[AmatT];
AmatTT=Prepend[AmatTT,datafc];
AmatTT=Prepend[AmatTT,datacn];

Bmat=avgBmatrix;
BmatT=Transpose[Bmat];
BmatT=Prepend[BmatT,Table["LCP mV",{i,1,21}]];
BmatT=Prepend[BmatT,{143800,72080,36120,18110,9075,4549,2280,1144,573.7,288.1,145.0,73.23,37.28,19.26,10.23,5.702,3.434,2.297,143800,0.290,0.290}];
BmatT=Prepend[BmatT,{0,3,6,9,12,15,18,21,24,27,30,33,36,39,42,45,48,51,0,"50 ohms on MC","Kay @ 101 dB"}];
BmatTT=Transpose[BmatT];
BmatTT=Prepend[BmatTT,datafc];
BmatTT=Prepend[BmatTT,datacn];

datamat=Table[0,{j,1,44},{i,1,chartchannels+3}];
datamat[[1;;3]]=AmatTT[[1;;3]];
Do[{
datamat[[k+2]]=AmatTT[[k+2]];
datamat[[k+2+21]]=BmatTT[[k+2]];
}
,{k,1,21}
];

datamatT=Transpose[datamat];
datamatT=Prepend[datamatT,"Values below in ADC millivolts, averaged from "<>ToString[sweeps]<>" sweeps ("<>IntegerString[Round2[chartduration]]<>" seconds of observation)"];
datamatT=Prepend[datamatT,"Temperatures are at DPS input terminals and do not include receiver noise figure"];
datamatT=Prepend[datamatT,"Noise source: HP 461 @ 68.2 MK through splitter with 0.26 dB measured insertion loss"];
datamatT=Prepend[datamatT,"Frequency sweep range "<>ToString[freqlo]<>" to "<>ToString[freqhi]<>" MHz"];
datamatT=Prepend[datamatT,"Step calibration performed on "<>chartstartF];
datamatT=Prepend[datamatT,"CSV file created "<>DateString[{"Day","-","MonthNameShort","-","Year"," ","Hour24",":","Minute"," Local"}]];
datamatT=Prepend[datamatT,csvfilename];
datamatT=Prepend[datamatT,"DPS Step Calibration Data"];

Export[csvfilename,datamatT]

DPS Color Scale Maker

ColorOffset=1550;
ColorGain=1.8;
ScaleWidth=4096;

StepCalInFile=Import["N:\\Documents\\0 RA\\Hardware\\DPS\\20130520 DPS Step Cal\\00 DPS Step Calibration 20130520.csv","Data"];

StepCalInFile[[All,2]]=(StepCalInFile[[All,2]]-ColorOffset)*ColorGain;

ListLogLinearPlot[StepCalInFile]
StepCalInFile//MatrixForm

iFn[x_]:=Interpolation[StepCalInFile][x];

Print["For a scale height of ",ScaleWidth," pixels:"]

Print["5 kK at ",iFn[5]]
Print["10 kK at ",iFn[10]]
Print["20 kK at ",iFn[30]]
Print["50 kK at ",iFn[50]]

Print["100 kK at ",iFn[100]]
Print["200 kK at ",iFn[200]]
Print["500 kK at ",iFn[500]]

Print["1 MK at ",iFn[1000]]
Print["2 MK at ",iFn[2000]]
Print["5 MK at ",iFn[5000]]

Print["10 MK at ",iFn[10000]]
Print["20 MK at ",iFn[20000]]
Print["50 MK at ",iFn[50000]]

Print["100 MK at ",iFn[100000]]
Print["200 MK at ",iFn[200000]]

ColorInFile=Import["C:\\Program Files\\Spectrograph\\AJ4CO-Rainbow.txt","Data"];
lo=1;
hi=Dimensions[ColorInFile][[1]];
ColorsIn ={
FractionalPart[FractionalPart[ColorInFile[[All,1]]/65536]*65536/256]*256,
IntegerPart[FractionalPart[ColorInFile[[All,1]]/65536]*65536/256],
IntegerPart[ColorInFile[[All,1]]/65536]
};
ColorsInT=Transpose[ColorsIn];
ColorsInT=ColorsInT/255;
xticks={
{310.`,Rotate["10 kK",-π/2],{0.01`,0},{}},
{434.`,Rotate["20 kK",-π/2],{0.01`,0},{}},
{539.`,Rotate["50 kK",-π/2],{0.01`,0},{}},
{712.`,Rotate["100 kK",-π/2],{0.01`,0},{}},
{917.`,Rotate["200 kK",-π/2],{0.01`,0},{}},
{1212.`,Rotate["500 kK",-π/2],{0.01`,0},{}},
{1467.`,Rotate["1.0 MK",-π/2],{0.01`,0},{}},
{1722.`,Rotate["2.0 MK",-π/2],{0.01`,0},{}},
{2010.`,Rotate["5.0 MK",-π/2],{0.01`,0},{}},
{2294.`,Rotate["10 MK",-π/2],{0.01`,0},{}},
{2558.`,Rotate["20 MK",-π/2],{0.01`,0},{}},
{2947.`,Rotate["50 MK",-π/2],{0.01`,0},{}},
{3212.`,Rotate["100 MK",-π/2],{0.01`,0},{}},
{3473.`,Rotate["200 MK",-π/2],{0.01`,0},{}},
{3932.`,Rotate["500 MK",-π/2],{0.01`,0},{}}
};

p=Rotate[Graphics[Raster[{ColorsInT}],AspectRatio->.07,Frame->True,FrameLabel->"XTFD Antenna Temperature",LabelStyle->Directive[FontFamily->"Arial",FontSize->12],FrameTicks->{xticks,None,None,None},Background->RGBColor[240/255,240/255,243/255],ImageSize->525],π/2]

Export["N:\\Documents\\0 RA\\Hardware\\DPS\\20130520 DPS Step Cal\\DPS Temp Scale Offset=1550 Gain=1.8.png",p]

ColorOffset=1550;
ColorGain=1.8;
ScaleWidth=4096;

StepCalInFile=Import["N:\\Documents\\0 RA\\Hardware\\DPS\\20130520 DPS Step Cal\\00 DPS Step Calibration 20130520.csv","Data"];


StepCalInFile//MatrixForm

StepCalInFile[[All,2]]=(StepCalInFile[[All,2]]-ColorOffset)*ColorGain;

StepCalInFile//MatrixForm

ListLogLinearPlot[StepCalInFile]
iFn[x_]:=Interpolation[StepCalInFile][x];

iFn[50]

DPS Correction Array Seed Generator - Takes short SPS file and saves continuum averages in CSV format

(********************************************************************************************)
(*                                                                                          *)
(*  Takes an SPS file snippet of the GB and calculates channel offsets necessary to         *)
(*  correct the observed GB to the theoretical GB model (see "The Galactic Background       *)
(*  in the Upper HF Band", Typinski, 2013).                                                 *)
(*                                                                                          *)
(*  Saves two CSV files, each an averaged continuum.                                        *)
(*                                                                                          *)
(*  File *_Correction_Array_Seed_Values.csv is the starting point for generating            *)
(*  the actual correction array using Excel (Excel is easier to work with to decide what    *)
(*  is bona-fide signal and what is noise or spurs and what is SWR ripple in the continua). *)
(*                                                                                          *)
(*  This routine uses the older method of parsing and displaying spectrogram data.          *)
(*  Decided to leave it this way as a ready reference.                                      *)
(*                                                                                          *)
(********************************************************************************************)

(* Set file name of input file                                  *)
filename="R:\\Observation Records\\AJ4CO_DPS\\1412\\141218\\141218065931 CA Snip.sps";

(* Set color file path and name for RCP & LCP spectrograms      *)
cpcolorfile="C:\\Program Files\\Spectrograph\\AJ4CO-Rainbow.txt";

(* Set color offset and color gain for RCP and LCP spectrograms *)
coloroffset=1975;
colorgain=2.5;

(* Create output filenames *)
outfile=StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1];
cCfilename=outfile<>"_Correction_Array_Seed_Values.csv";


(* Open one file for input                                      *)
infile=OpenRead[filename,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];
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=(freqhi-freqlo)/chartchannels*10^3;
channeldwelltime=Round2[chartduration/(chartrecs*chartchannels)*10^6];
sweeps=ToExpression[StringTake[chartnote,{ StringPosition[chartnote,"SWEEPS"][[1,2]]+1 ,StringPosition[chartnote,"ÿLOWF"][[1,1]]-1}]];

(* Create temp table and read in spectro data.                  *)
(* 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)                *)
indata=Table[0,{i,1,sweeps},{j,1,chartchannels*2+1}];
Do[{
indata[[i,j]]=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
}
,{i,1,sweeps}
,{j,1,chartchannels*2+1}
];
Close[infile];

(* Verify number of records present in data file by checking    *)
(* count of EOR markers.                                        *)
(* The EOR marker is two consecutive 0xFE characters, ASCII "þ" *)
(* or 65278 decomal when read as a 2-byte unsigned integer.     *)
chartrecs=Count[indata,65278,2];

(* Calculate chart sweep rate                                   *)
chartsweeprate=chartrecs/chartduration;

(* Find min and max of raw data over all data points            *)
chartmaxF=Max[indata[[1;;chartrecs,1;;chartchannels*2]]];
chartminF=Min[indata[[1;;chartrecs,1;;chartchannels*2]]];

(* Print formatted file header data to screen                   *)
Print["DPS Correction Array Seed Generator"]
Print["Analysis run date: "<>DateString[{"Day","-","MonthNameShort","-","Year"," ","Hour24",":","Minute"," Local"}]]
Print[]
Print["         Local name: ",chartlocalname]
Print["         Input file: ", filename]
Print["        RSS version: ",chartversionF]
Print["        Chart start: ",chartstartF]
Print["          Chart end: ",chartendF]
Print[]
Print["         Local name: ",chartlocalname]
Print["     Chart location: ",chartlocation]
Print["           Latitude: ",ToString[PaddedForm[chartlat,{8,6},NumberPadding->{"0","0"},NumberSigns->{"(S) -","(N) +"}]]<>" degrees"]
Print["          Longitude: ",ToString[PaddedForm[chartlong,{8,6},NumberPadding->{"0","0"},NumberSigns->{"(W) -","(E) +"}]]<>" degrees"]
Print["          Time zone: UTC+",chartTZ]
Print["           Observer: ",chartauthor]
Print[]
Print["        Data Source: ",chartsource]
Print["      Data duration: ",ToString[PaddedForm[chartduration,{10,3},NumberPadding->{"","0"}]]," s"]
Print["        Sweep Start: ",ToString[PaddedForm[freqhi,{10,6},NumberPadding->{"","0"}]]," MHz"]
Print["          Sweep End: ",ToString[PaddedForm[freqlo,{10,6},NumberPadding->{"","0"}]]," MHz"]
Print["             Sweeps: ",chartrecs," ",sweeps]
Print["         Sweep rate: ",ToString[PaddedForm[chartsweeprate,{6,3},NumberPadding->{"","0"}]]," Hz"]
Print["           Channels: ",chartchannels]
Print[" Channel separation: ",ToString[PaddedForm[channelsep*1.,{7,4},NumberPadding->{"","0"}]]," kHz"]
Print[" Channel dwell time: ",ToString[channeldwelltime]<>" μs"]
Print[]
Print["           Data Max: ",ToString[chartmaxF]<>" mV"]
Print["           Data Min: ",ToString[chartminF]<>" mV"]
Print[]
Print["  Chart note length: ",chartnotelength]
Print["       Chart note 1: ",chartnote]

(* Create 12-bit data from input file array; this strips out    *)
(* the EOR markers from the array.                              *)
(* DPS SPS files are already 12-bit, so no multiplication is    *)
(* necessary.                                                   *)
indata12bit=indata[[1;;chartrecs,1;;chartchannels*2]];

(* Create DPS Side A and Side B arrays.                         *)
(* Normal operation: Side A = RCP, Side B = LCP                 *)
(* Extract RCP and LCP data from interleaved raw data           *)
(* Raw data record format:                                      *)
(* Side B1, Side A1, Side B2, Side A2 ... Side Bn, Side An, EOR *)
dataA=Table[0.0,{i,1,chartrecs},{j,1,chartchannels}];
dataB=Table[0.0,{i,1,chartrecs},{j,1,chartchannels}];
Do[
Do[
{dataA[[i,j]]=indata12bit[[i,j*2]],
dataB[[i,j]]=indata12bit[[i,j*2-1]]}
,{j,1,chartchannels}
],{i,1,chartrecs}
];

(*  RSS and DPS data have channel 200 rolled over above channel 1. *)
(*  These two lines put channel 200 data back where it belongs.    *)
dataA=RotateLeft[dataA,{0,1}];
dataB=RotateLeft[dataB,{0,1}];

(* Apply color offset and color gain to raw data, just like RSS *) 
dataAadj=Round2[(dataA-coloroffset)*colorgain];
dataBadj=Round2[(dataB-coloroffset)*colorgain];

(* Place-holder for a data shrinkage routine                    *)
(* Not really needed for small files (a minute or two of data)  *)
dataAshrunk=dataAadj;
dataBshrunk=dataBadj;
chartshrunkrecs=chartrecs;

(* Open and read color file data for circular pol spectrograms  *)
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;

(* Create raster arrays of RGB triplets of raw spectro data     *)
dataAshrunkraster=dataAshrunk;
dataBshrunkraster=dataBshrunk;
(* Replace raw data values with color triplets based on color   *)
(* file data                                                    *)
(* Colorize DPS Side A shrunk data                              *)  
Do[{
If[
dataAshrunkraster[[i,j]]≤colorhi,
{If[
dataAshrunkraster[[i,j]]≥colorlo,
{dataAshrunkraster[[i,j]]=ColorsInT[[dataAshrunkraster[[i,j]]]]},
{dataAshrunkraster[[i,j]]=ColorsInT[[colorlo]]}
]
},
{dataAshrunkraster[[i,j]]=ColorsInT[[colorhi]]}
]
}
,{i,1,chartshrunkrecs}
,{j,1,chartchannels}
]
(* Colorize DPS Side B shrunk data                              *)
Do[{
If[
dataBshrunkraster[[i,j]]≤colorhi,
{If[
dataBshrunkraster[[i,j]]≥colorlo,
{dataBshrunkraster[[i,j]]=ColorsInT[[dataBshrunkraster[[i,j]]]]},
{dataBshrunkraster[[i,j]]=ColorsInT[[colorlo]]}
]
},
{dataBshrunkraster[[i,j]]=ColorsInT[[colorhi]]}
]
}
,{i,1,chartshrunkrecs}
,{j,1,chartchannels}
]
(* Transpose arrays to match Mathematica's default for raster   *)
(* images.                                                      *)
dataAshrunkrasterT=Reverse[Transpose[dataAshrunkraster],1];
dataBshrunkrasterT=Reverse[Transpose[dataBshrunkraster],1];

(* Print text info and spectrograms                             *)
Print[]
Print[]
Print["Raw Data Analysis"]
Print[]
Print["         Color File: "<>cpcolorfile]
Print["       Color Offset: ",IntegerString[coloroffset]]
Print["         Color Gain: ",ToString[PaddedForm[colorgain,{2,1},NumberPadding->{"","0"}]]]
Print["       Decimated to: ",IntegerString[chartshrunkrecs]," sweeps"]
Print[]

Print["       Polarization: RCP - raw data"]
ImageResize[Image[Raster[dataAshrunkrasterT, Automatic]],{733,300},Resampling-> "Bilinear"]
Print["       Polarization: LCP - raw data"]
Image[Raster[dataBshrunkrasterT, Automatic]]


(* Calculate time averaged continua for Side A & Side B data    *)
dataAM=N[Mean[dataA]];
dataBM=N[Mean[dataB]];

(* Print chart label and continua plot                          *)
Print[]
Print["      Continuum display (mean values of ",ToString[PaddedForm[chartduration,{10,3},NumberPadding->{"","0"}]]," seconds of data)"]
ListLinePlot[{dataAM,dataBM},PlotRange->{1500,3000},PlotStyle->{Darker[Green],Red},AspectRatio->0.6]

(* Create channel frequency list                                *)
dataFM=Table[ToString[PaddedForm[N[freqhi-i*channelsep*10^-3],{6,4},NumberPadding->{"","0"}]],{i,1,chartchannels}];

(* Export processed data to output files                        *)
dataCM=Table[{i,0,0,0},{i,1,chartchannels}];
dataCM[[All,2]]=dataFM;
dataCM[[All,3]]=Round2[dataAM];
dataCM[[All,4]]=Round2[dataBM];
dataCM=Prepend[dataCM,{"Channel","f_c (MHz)","Side A (mV)","Side B (mV)"}];
dataCM=Prepend[dataCM,"Values below averaged from "<>ToString[sweeps]<>" sweeps ("<>ToString[PaddedForm[chartduration,{10,3},NumberPadding->{"","0"}]]<>" seconds of observation)"];
dataCM=Prepend[dataCM,"Frequency sweep range "<>ToString[freqlo]<>" to "<>ToString[freqhi]<>" MHz"];
dataCM=Prepend[dataCM,"Continuum data from "<>filename];
dataCM=Prepend[dataCM,"CSV file created "<>DateString[{"Day","-","MonthNameShort","-","Year"," ","Hour24",":","Minute"," Local"}]];
dataCM=Prepend[dataCM,cCfilename];
dataCM=Prepend[dataCM,"DPS Correction Array Seed Values"];
Export[cCfilename,dataCM]

DPS File Correction Array Applyer  - 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="R:\\Observation Records\\AJ4CO_DPS\\1501\\150116\\150116054500.sps";

(* Define correction array to be applied *)
carrfilename="R:\\Observation Records\\AJ4CO_DPS\\1412\\141218\\DPS Correction Array 2014 12 18 B.csv";



infile=OpenRead[filename,BinaryFormat->True];
outfile=OpenWrite[StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>" corrected using CA 2014 12 18 B.sps",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
];

(* 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]];
(* Write data to output file                                 *)
Do[
{
BinaryWrite[outfile,indataB[[i]],"UnsignedInteger16",ByteOrdering->+1];
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];
}
,{j,1,sweeps}
];
Close[infile];
Close[outfile]
Print["Done."]

DPS Data Analyzer - 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.                     *)
(*                                                                                          *)
(********************************************************************************************)

(* Set file name of input file                                  *)
filename="R:\\Observation Records\\AJ4CO_DPS-XTFD\\1307\\130705\\130706015208.sps";

(* Set file name of correction array                            *)
cafilename="R:\\Observation Records\\AJ4CO_DPS-XTFD\\1305\\130522\\130522_Correction_Array.csv";

(* Set 20 MHz GB cal point in mV and detector rate in mV/dB *)
calpoint=1850;
detrate=45;

(* Set color offset and color gain for RCP and LCP spectrograms *)
coloroffset=1550;
colorgain=1.8;

(* Set Xpol spectrogram color gain and color offset             *)
xcolorgain=3;
xcoloroffset=-2048;

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

(* Set color file path and name for RCP & LCP spectrograms      *)
cpcolorfile="C:\\Program Files\\Spectrograph\\AJ4CO-Rainbow.txt";

(* Set color file path and name for X-polarization spectrogram  *)
xpcolorfile="C:\\Program Files\\Spectrograph\\CrossPolarization.txt";

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

(* Set strip chart channel number *)
stripchn=161;

(* 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]];

(* Create output filenames *)
rptfilename=StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>"_spectrogram_inspection.pdf";
spgmfilename=StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>"_spectrogram_DCP_raw.png";
spgmAfilename=StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>"_spectrogram_RCP_raw.png";
spgmBfilename=StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>"_spectrogram_LCP_raw.png";
spgmXfilename=StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>"_spectrogram_X-pol_raw.png";
contfilename=StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>"_continua_raw.png";
stripfilename=StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>"_strip_chart_channel_"<>ToString[stripchn]<>"_raw.png";
spgmfilenamecorr=StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>"_spectrogram_DCP_corrected.png";
spgmAfilenamecorr=StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>"_spectrogram_RCP_corrected.png";
spgmBfilenamecorr=StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>"_spectrogram_LCP_corrected.png";
spgmXfilenamecorr=StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>"_spectrogram_X-pol_corrected.png";
contfilenamecorr=StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>"_continua_corrected.png";
stripfilenamecorr=StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>"_strip_chart_channel_"<>ToString[stripchn]<>"_corrected.png";


(* 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;
(* Graphics[Raster[{ColorsInT}],AspectRatio->.2] *)

(* Repeat above for Xpol color file                             *)
xColorInFile=Import[xpcolorfile,"Data"];
xcolorlo=1;
xcolorhi=Dimensions[xColorInFile][[1]];
xColorsIn ={
FractionalPart[FractionalPart[xColorInFile[[All,1]]/65536]*65536/256]*256,
IntegerPart[FractionalPart[xColorInFile[[All,1]]/65536]*65536/256],
IntegerPart[xColorInFile[[All,1]]/65536]
};
xColorsInT=Transpose[xColorsIn];
ColorsInTx=xColorsInT/255;
(* Graphics[Raster[{ColorsInTx}],AspectRatio->.2] *)

(* Open one SPS file for input                                  *)
infile=OpenRead[filename,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 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 >"];

(* Create temp table and read in spectro data.                  *)
(* 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)                *)
dataAdecimated=Table[0,{i,1,chartwidth},{j,1,chartchannels}];
dataBdecimated=Table[0,{i,1,chartwidth},{j,1,chartchannels}];
(* Create input record accumulators to sum & average input data *)
indataArec=Table[0,{i,1,chartchannels}];
indataBrec=Table[0,{i,1,chartchannels}];
shrinkratio=sweeps/chartwidth;
(* Initialize raw data min-max accumulators                     *)
chartAmaxR=0;
chartAminR=2^nbits-1;
chartBmaxR=0;
chartBminR=2^nbits-1;
(* Read data for all ouput chart records                        *)
Do[
{
(* Write averaged input records to one chart record          *)
Do[
{
(* Read input data record                                  *)
Do[
{
(* Read Side B data (side B comes first in the SPS data) *)
tempB=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
chartBmaxR=Max[chartBmaxR,tempB];
chartBminR=Min[chartBminR,tempB];
indataBrec[[i]]=indataBrec[[i]]+tempB;
(* Read Side A data                                      *)
tempA=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
chartAmaxR=Max[chartAmaxR,tempA];
chartAminR=Min[chartAminR,tempA];
indataArec[[i]]=indataArec[[i]]+tempA;
}
(* cycle over 2 x number of channels *)
,{i,1,chartchannels}
];
(* Read EOR marker                                         *)
dummy=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
}
(* cycle over number of input sweeps to be averaged         *)
,{j,Round2[1+(z-1)*(shrinkratio)],Round2[z*shrinkratio]} 
];
(*  RSS and DPS data have channel 200 rolled over above channel 1. *)
(*  These two lines put channel 200 data back where it belongs.    *)
indataArec=RotateLeft[indataArec];
indataBrec=RotateLeft[indataBrec];
(* Average summed input records and write to decimated array *)
(* See "Subsampling Area Averaging Method.xls" for details   *)
dataAdecimated[[z]]=Round2[(indataArec/(Round2[z*shrinkratio]-Round2[1+(z-1)*shrinkratio]+1))];
dataBdecimated[[z]]=Round2[(indataBrec/(Round2[z*shrinkratio]-Round2[1+(z-1)*shrinkratio]+1))];
(* Zero out input record accumulators before next round of averaging *)
indataArec=indataArec*0;
indataBrec=indataBrec*0;
}
(* cycle over output chart width in pixels, sweeps, records, whatever *)
,{z,1,chartwidth}
];
Close[infile];

(* Apply correction arrays to decimated data                    *)
dataAdecimatedcorr=dataAdecimated;
dataBdecimatedcorr=dataBdecimated;
Do[dataAdecimatedcorr[[i]]=dataAdecimated[[i]]+corrarray[[1]],{i,1,chartwidth}];
Do[dataBdecimatedcorr[[i]]=dataBdecimated[[i]]+corrarray[[2]],{i,1,chartwidth}];

(* Find min and max of decimated data over all data points      *)
chartAmaxD=Round2[Max[dataAdecimated]];
chartAminD=Round2[Min[dataAdecimated]];
chartBmaxD=Round2[Max[dataBdecimated]];
chartBminD=Round2[Min[dataBdecimated]];
chartAmaxDcorr=Round2[Max[dataAdecimatedcorr]];
chartAminDcorr=Round2[Min[dataAdecimatedcorr]];
chartBmaxDcorr=Round2[Max[dataBdecimatedcorr]];
chartBminDcorr=Round2[Min[dataBdecimatedcorr]];

(* Apply color offset and color gain to raw data, just like RSS *) 
dataAadj=Round2[(dataAdecimated-coloroffset)*colorgain];
dataBadj=Round2[(dataBdecimated-coloroffset)*colorgain];
dataAadjcorr=Round2[(dataAdecimatedcorr-coloroffset)*colorgain];
dataBadjcorr=Round2[(dataBdecimatedcorr-coloroffset)*colorgain];

(* Coerce adjusted data that lays outside n-bit range          *)
datalo=1;
datahi=2^nbits-1;  
(* Coerce Side A data                                          *)
Do[{
If[
dataAadj[[i,j]]≤datahi,
{If[
dataAadj[[i,j]]≥datalo,
{dataAadj[[i,j]]=dataAadj[[i,j]]},
{dataAadj[[i,j]]=datalo}
]
},
{dataAadj[[i,j]]=datahi}
]
}
,{i,1,chartwidth}
,{j,1,chartchannels}
]
(* Coerce Side B data                                          *)
Do[{
If[
dataBadj[[i,j]]≤datahi,
{If[
dataBadj[[i,j]]≥datalo,
{dataBadj[[i,j]]=dataBadj[[i,j]]},
{dataBadj[[i,j]]=datalo}
]
},
{dataBadj[[i,j]]=datahi}
]
}
,{i,1,chartwidth}
,{j,1,chartchannels}
]

(* Coerce corrected adjusted data that lays outside n-bit range *)
datalo=1;
datahi=2^nbits-1;  
(* Coerce Side A data                                          *)
Do[{
If[
dataAadjcorr[[i,j]]≤datahi,
{If[
dataAadjcorr[[i,j]]≥datalo,
{dataAadjcorr[[i,j]]=dataAadjcorr[[i,j]]},
{dataAadjcorr[[i,j]]=datalo}
]
},
{dataAadjcorr[[i,j]]=datahi}
]
}
,{i,1,chartwidth}
,{j,1,chartchannels}
]
(* Coerce Side B data                                          *)
Do[{
If[
dataBadjcorr[[i,j]]≤datahi,
{If[
dataBadjcorr[[i,j]]≥datalo,
{dataBadjcorr[[i,j]]=dataBadjcorr[[i,j]]},
{dataBadjcorr[[i,j]]=datalo}
]
},
{dataBadjcorr[[i,j]]=datahi}
]
}
,{i,1,chartwidth}
,{j,1,chartchannels}
]

(* Create X-pol data array                                     *)
dataXdecimated=dataAdecimated;
dataXdecimated=dataAdecimated-dataBdecimated;
xchartmax=Round2[Max[dataXdecimated]];
xchartmin=Round2[Min[dataXdecimated]];
dataXadj=Round2[dataXdecimated*xcolorgain-xcoloroffset];
(* Coerce adjusted data that lays outside n-bit range          *)
datalo=1;
datahi=2^nbits-1;  
(* Coerce X-pol data                                           *)
Do[{
If[
dataXadj[[i,j]]≤datahi,
{If[
dataXadj[[i,j]]≥datalo,
{dataXadj[[i,j]]=dataXadj[[i,j]]},
{dataXadj[[i,j]]=datalo}
]
},
{dataXadj[[i,j]]=datahi}
]
}
,{i,1,chartwidth}
,{j,1,chartchannels}
]

(* Create corrected X-pol data array                           *)
dataXdecimatedcorr=dataAdecimatedcorr;
dataXdecimatedcorr=dataAdecimatedcorr-dataBdecimatedcorr;
xchartmaxcorr=Round2[Max[dataXdecimatedcorr]];
xchartmincorr=Round2[Min[dataXdecimatedcorr]];
dataXadjcorr=Round2[dataXdecimatedcorr*xcolorgain-xcoloroffset];
(* Coerce adjusted data that lays outside n-bit range          *)
datalo=1;
datahi=2^nbits-1;  
(* Coerce X-pol data                                           *)
Do[{
If[
dataXadjcorr[[i,j]]≤datahi,
{If[
dataXadjcorr[[i,j]]≥datalo,
{dataXadjcorr[[i,j]]=dataXadjcorr[[i,j]]},
{dataXadjcorr[[i,j]]=datalo}
]
},
{dataXadjcorr[[i,j]]=datahi}
]
}
,{i,1,chartwidth}
,{j,1,chartchannels}
]

(* Calculate time averaged continua for Side A & Side B data    *)
dataAM=N[Mean[dataAdecimated]];
dataBM=N[Mean[dataBdecimated]];
dataAMcorr=N[Mean[dataAdecimatedcorr]];
dataBMcorr=N[Mean[dataBdecimatedcorr]];

(* Define galactic background model, output in kK               *)
GB[f_]:=1.07*10^5*f^-2.56;
(* Define function to convert GB kK to ADC mV                   *)
GBmV[f_]:=10 Log[10,GB[f]/GB[20]]*detrate+calpoint;
(* Create table of GB model ADC values in mV                    *)
dataGB=Table[GBmV[f],{f,33-0.08,17,-0.08}];

(* Print formatted file header data to screen                   *)
output=Table["",{i,1,150}];
output[[1]]=Style["DPS Data Analysis",Larger,Darker[Darker[Blue]],Bold];
output[[2]]="Run date: "<>DateString[{"Day","-","MonthNameShort","-","Year"," ","Hour24",":","Minute"," Local"}];

output[[10]]=Style["File Information – Raw Data",Bold];

output[[12]]="         Local name:  "<>chartlocalname;
output[[13]]="         Input file:  "<> filename;
output[[14]]="        RSS version:  "<>chartversionF;
output[[15]]="        Chart start:  "<>chartstartF;
output[[16]]="          Chart end:  "<>chartendF;

output[[18]]="         Local name:  "<>chartlocalname;
output[[19]]="     Chart location:  "<>chartlocation;
output[[20]]="           Latitude:  "<>ToString[PaddedForm[chartlat,{8,6},NumberPadding->{"0","0"},NumberSigns->{"(S) -","(N) +"}]]<>" degrees";
output[[21]]="          Longitude:  "<>ToString[PaddedForm[chartlong,{8,6},NumberPadding->{"0","0"},NumberSigns->{"(W) -","(E) +"}]]<>" degrees";
output[[22]]="          Time zone:  UTC+"<>ToString[chartTZ];
output[[23]]="           Observer:  "<>chartauthor;

output[[25]]="        Data source:  "<>chartsource;
output[[26]]="      Data duration:  "<>ToString[PaddedForm[chartduration,{10,3},NumberPadding->{"","0"}]]<>" s";
output[[27]]="        Sweep start:  "<>ToString[PaddedForm[freqhi,{10,3},NumberPadding->{"","0"}]]<>" MHz";
output[[28]]="          Sweep end:  "<>ToString[PaddedForm[freqlo,{10,3},NumberPadding->{"","0"}]]<>" MHz";
output[[29]]="             Sweeps:  "<>ToString[sweeps];
output[[30]]="         Sweep rate:  "<>ToString[PaddedForm[chartsweeprate,{6,3},NumberPadding->{"","0"}]]<>" Hz";
output[[31]]="           Channels:  "<>ToString[chartchannels];
output[[32]]=" Channel separation:  "<>ToString[PaddedForm[channelsep,{6,3},NumberPadding->{"","0"}]]<>" kHz";
output[[33]]=" Channel dwell time:  "<>ToString[channeldwelltime]<>" μs";

output[[35]]="RCP ADC range (raw):  "<>ToString[chartAminR]<>" mV min  to  "<>ToString[chartAmaxR]<>" mV max";
output[[36]]="LCP ADC range (raw):  "<>ToString[chartBminR]<>" mV min  to  "<>ToString[chartBmaxR]<>" mV max";

output[[38]]="   Data note length:  "<>ToString[chartnotelength];
output[[39]]="          Data note:  "<>chartnote;

output[[41]]="   Correction array:  "<>cafilename;

output[[43]]=Style["Spectrograms – Decimated                ",Bold];

output[[45]]="         Color file:  "<>cpcolorfile;
output[[46]]="       Color offset:  "<>ToString[coloroffset];
output[[47]]="         Color gain:  "<>ToString[PaddedForm[colorgain,{2,1},NumberPadding->{"","0"}]];
output[[48]]="               Data:  "<>ToString[sweeps]<>" sweeps";
output[[49]]="       Decimated to:  "<>ToString[chartwidth]<>" sweeps  (decimation factor "<>ToString[PaddedForm[N[shrinkratio],{10,3},NumberPadding->{"","0"}]]<>" )";
output[[50]]="  Decimation method:  "<>"Mean (arithmetic average)";
output[[51]]="RCP ADC range (decimated):  "<>ToString[chartAminD]<>" mV min  to  "<>ToString[chartAmaxD]<>" mV max";
output[[52]]="LCP ADC range (decimated):  "<>ToString[chartBminD]<>" mV min  to  "<>ToString[chartBmaxD]<>" mV max";

plot1=ArrayPlot[Transpose[dataAadj],
PixelConstrained->False,
ImageSize->{1000,300},
ColorFunctionScaling->False,
ColorFunction->Function[{z},RGBColor[ColorsInT[[z]]]],
PlotRangePadding->{{10,10},{10,10}},
Background->None,
Frame->True,
FrameTicks->True,
RotateLabel-> True,
FrameStyle->{{Directive[Thin,FontFamily->Helvetica,11],Directive[Thin,FontFamily->Helvetica,11]},{Directive[Thin,FontFamily->Helvetica,11],Directive[Thin,FontFamily->Helvetica,11]}},FrameLabel->{Style["RCP   Channel",FontSize->11,FontFamily->Helvetica],Style["",FontSize->11,FontFamily->Helvetica],"RCP   Frequency (MHz)",""}
];

plot2=ArrayPlot[Transpose[dataBadj],
PixelConstrained->False,
ImageSize->{1000,300},
ColorFunctionScaling->False,
ColorFunction->Function[{z},RGBColor[ColorsInT[[z]]]],
PlotRangePadding->{{10,10},{10,10}},
Background->None,
Frame->True,
FrameTicks->True,
RotateLabel-> True,
FrameStyle->{{Directive[Thin,FontFamily->Helvetica,11],Directive[Thin,FontFamily->Helvetica,11]},{Directive[Thin,FontFamily->Helvetica,11],Directive[Thin,FontFamily->Helvetica,11]}},FrameLabel->{Style["LCP   Channel",FontSize->11,FontFamily->Helvetica],Style["",FontSize->11,FontFamily->Helvetica],"LCP   Frequency (MHz)",""}
];
back=Image[Raster[Table[{240/255,240/255,243/255},{i,1,546},{j,1,1000}]],ImageSize->{1000,546}];
mid=ImageCompose[back,plot1,{500,396}];
dual=ImageCompose[mid,plot2,{500,150}];
output[[54]]=dual;

output[[56]]=Style["Spectrograms – Decimated – Correction Arrays Applied",Bold];

output[[58]]="         Color file:  "<>cpcolorfile;
output[[59]]="       Color offset:  "<>ToString[coloroffset];
output[[60]]="         Color gain:  "<>ToString[PaddedForm[colorgain,{2,1},NumberPadding->{"","0"}]];
output[[61]]="               Data:  "<>ToString[sweeps]<>" sweeps";
output[[62]]="       Decimated to:  "<>ToString[chartwidth]<>" sweeps  (decimation factor "<>ToString[PaddedForm[N[shrinkratio],{10,3},NumberPadding->{"","0"}]]<>" )";
output[[63]]="  Decimation method:  "<>"Mean (arithmetic average)";
output[[64]]="RCP ADC range (decimated):  "<>ToString[chartAminDcorr]<>" mV min  to  "<>ToString[chartAmaxDcorr]<>" mV max";
output[[65]]="LCP ADC range (decimated):  "<>ToString[chartBminDcorr]<>" mV min  to  "<>ToString[chartBmaxDcorr]<>" mV max";

plot6=ArrayPlot[Transpose[dataAadjcorr],
PixelConstrained->False,
ImageSize->{1000,300},
ColorFunctionScaling->False,
ColorFunction->Function[{z},RGBColor[ColorsInT[[z]]]],
PlotRangePadding->{{10,10},{10,10}},
Background->None,
Frame->True,
FrameTicks->True,
RotateLabel-> True,
FrameStyle->{{Directive[Thin,FontFamily->Helvetica,11],Directive[Thin,FontFamily->Helvetica,11]},{Directive[Thin,FontFamily->Helvetica,11],Directive[Thin,FontFamily->Helvetica,11]}},FrameLabel->{Style["RCP   Channel",FontSize->11,FontFamily->Helvetica],Style["",FontSize->11,FontFamily->Helvetica],"RCP   Frequency (MHz)",""}
];

plot7=ArrayPlot[Transpose[dataBadjcorr],
PixelConstrained->False,
ImageSize->{1000,300},
ColorFunctionScaling->False,
ColorFunction->Function[{z},RGBColor[ColorsInT[[z]]]],
PlotRangePadding->{{10,10},{10,10}},
Background->None,
Frame->True,
FrameTicks->True,
RotateLabel-> True,
FrameStyle->{{Directive[Thin,FontFamily->Helvetica,11],Directive[Thin,FontFamily->Helvetica,11]},{Directive[Thin,FontFamily->Helvetica,11],Directive[Thin,FontFamily->Helvetica,11]}},FrameLabel->{Style["LCP   Channel",FontSize->11,FontFamily->Helvetica],Style["",FontSize->11,FontFamily->Helvetica],"LCP   Frequency (MHz)",""}
];
back=Image[Raster[Table[{240/255,240/255,243/255},{i,1,546},{j,1,1000}]],ImageSize->{1000,546}];
mid=ImageCompose[back,plot6,{500,396}];
dual=ImageCompose[mid,plot7,{500,150}];
output[[67]]=dual;

output[[69]]=Style["Continua – Decimated",Bold];
(* Create a continuum plot                                      *)
charttitle=Column[{"",
Style["RCP & LCP Continua – average of "<>ToString[chartwidth]<>" sweeps"<>" decimated from "<>ToString[sweeps]<>" sweeps  ("<>ToString[Round2[chartduration]]<>" seconds in original data)",Plain,FontSize->12,FontFamily->Helvetica],
Row[{
Style["RCP = Green",Plain,FontSize->12,FontFamily->Helvetica,RGBColor[0,150/255,100/255]],
Style["LCP = RED",Plain,FontSize->12,FontFamily->Helvetica,RGBColor[190/255,39/255,39/255]],
Style["GB Model = Blue (for slope ONLY)",Plain,FontSize->12,FontFamily->Helvetica,RGBColor[0,0,160/255]]
},"        "]
},Center];
plot3=ListLinePlot[{dataBM,dataAM,dataGB},
DataRange->{1,chartchannels},
PlotRange->{1500,2100},
PlotStyle->{{Thick,RGBColor[190/255,39/255,39/255]},{Thick,RGBColor[0,150/255,100/255]},{Thick,RGBColor[0,0,160/255]}},
AspectRatio->200/chartwidth,
PlotRangePadding->{{3,3},{3,3}},
PlotMarkers->None,
Axes->None,
GridLines->Automatic,
Background->RGBColor[240/255,240/255,243/255],
Frame->True,
RotateLabel-> True,
FrameStyle->{{Directive[Thin,FontFamily->Helvetica,11],Thin},{Directive[Thin,FontFamily->Helvetica,11],Thin}},FrameLabel->{Style["Channel Number",FontSize->11,FontFamily->Helvetica],Style["Average ADC Value (mV)",FontSize->11,FontFamily->Helvetica],charttitle,""}
];
back=Image[Raster[Table[{240/255,240/255,243/255},{i,1,300},{j,1,1000}]],ImageSize->{1000,300}];
fore=Rasterize[plot3,ImageSize->chartwidth+110];
mid=ImageCompose[back,fore,{485,150}];
output[[70]]=mid;

output[[72]]=Style["Continua – Decimated – Correction Array Applied",Bold];
(* Create a continuum plot                                      *)
charttitle=Column[{"",
Style["RCP & LCP Continua – average of "<>ToString[chartwidth]<>" sweeps"<>" decimated from "<>ToString[sweeps]<>" sweeps  ("<>ToString[Round2[chartduration]]<>" seconds in original data)",Plain,FontSize->12,FontFamily->Helvetica],
Row[{
Style["RCP = Green",Plain,FontSize->12,FontFamily->Helvetica,RGBColor[0,150/255,100/255]],
Style["LCP = RED",Plain,FontSize->12,FontFamily->Helvetica,RGBColor[190/255,39/255,39/255]],
Style["GB Model = Blue (for slope ONLY)",Plain,FontSize->12,FontFamily->Helvetica,RGBColor[0,0,160/255]]
},"        "]
},Center];
plot8=ListLinePlot[{dataBMcorr,dataAMcorr,dataGB},
DataRange->{1,chartchannels},
PlotRange->{1500,2100},
PlotStyle->{{Thick,RGBColor[190/255,39/255,39/255]},{Thick,RGBColor[0,150/255,100/255]},{Thick,RGBColor[0,0,160/255]}},
AspectRatio->200/chartwidth,
PlotRangePadding->{{3,3},{3,3}},
PlotMarkers->None,
Axes->None,
GridLines->Automatic,
Background->RGBColor[240/255,240/255,243/255],
Frame->True,
RotateLabel-> True,
FrameStyle->{{Directive[Thin,FontFamily->Helvetica,11],Thin},{Directive[Thin,FontFamily->Helvetica,11],Thin}},FrameLabel->{Style["Channel Number",FontSize->11,FontFamily->Helvetica],Style["Average ADC Value (mV)",FontSize->11,FontFamily->Helvetica],charttitle,""}
];
back=Image[Raster[Table[{240/255,240/255,243/255},{i,1,300},{j,1,1000}]],ImageSize->{1000,300}];
fore=Rasterize[plot8,ImageSize->chartwidth+110];
mid=ImageCompose[back,fore,{485,150}];
output[[73]]=mid;

output[[75]]=Style["Strip Chart – Decimated",Bold];
(* Create a continuum plot                                      *)
charttitle=Column[{"",
Style["RCP & LCP Strip Chart – "<>ToString[chartwidth]<>" sweeps"<>" decimated from "<>ToString[sweeps]<>" sweeps – Channel "<>ToString[stripchn],Plain,FontSize->12,FontFamily->Helvetica],
Row[{
Style["RCP = Green",Plain,FontSize->12,FontFamily->Helvetica,RGBColor[0,150/255,100/255]],
Style["LCP = RED",Plain,FontSize->12,FontFamily->Helvetica,RGBColor[190/255,39/255,39/255]],
Style["GB Model = Blue (NOT YET MEANINGFUL)",Plain,FontSize->12,FontFamily->Helvetica,RGBColor[0,0,160/255]]
},"        "]
},Center];
plot4=ListLinePlot[{dataBdecimated[[All,stripchn]],dataAdecimated[[All,stripchn]],Table[dataGB[[stripchn]],{i,1,chartwidth}]},
DataRange->{1,chartwidth},
PlotRange->{1500,2500},
PlotStyle->{{Thick,RGBColor[190/255,39/255,39/255]},{Thick,RGBColor[0,150/255,100/255]},{Thick,RGBColor[0,0,160/255]}},
ImageSize->{chartwidth,400},
AspectRatio->200/chartwidth,
PlotRangePadding->{{3,3},{3,3}},
PlotMarkers->None,
Axes->None,
GridLines->Automatic,
Background->RGBColor[240/255,240/255,243/255],
Frame->True,
RotateLabel-> True,
FrameStyle->{{Directive[Thin,FontFamily->Helvetica,11],Thin},{Directive[Thin,FontFamily->Helvetica,11],Thin}},FrameLabel->{Style["Sweep Number (Decimated Data)",FontSize->11,FontFamily->Helvetica],Style["ADC Value (mV)",FontSize->11,FontFamily->Helvetica],charttitle,""}
];
fore=Rasterize[plot4,ImageSize->chartwidth+110];
mid=ImageCompose[back,fore,{485,150}];
output[[76]]=mid;

output[[78]]=Style["Strip Chart – Decimated – Correction Array Applied",Bold];
(* Create a continuum plot                                      *)
charttitle=Column[{"",
Style["RCP & LCP Strip Chart – "<>ToString[chartwidth]<>" sweeps"<>" decimated from "<>ToString[sweeps]<>" sweeps – Channel "<>ToString[stripchn],Plain,FontSize->12,FontFamily->Helvetica],
Row[{
Style["RCP = Green",Plain,FontSize->12,FontFamily->Helvetica,RGBColor[0,150/255,100/255]],
Style["LCP = RED",Plain,FontSize->12,FontFamily->Helvetica,RGBColor[190/255,39/255,39/255]],
Style["GB Model = Blue (NOT YET MEANINGFUL)",Plain,FontSize->12,FontFamily->Helvetica,RGBColor[0,0,160/255]]
},"        "]
},Center];
plot9=ListLinePlot[{dataBdecimatedcorr[[All,stripchn]],dataAdecimatedcorr[[All,stripchn]],Table[dataGB[[stripchn]],{i,1,chartwidth}]},
DataRange->{1,chartwidth},
PlotRange->{1500,2500},
PlotStyle->{{Thick,RGBColor[190/255,39/255,39/255]},{Thick,RGBColor[0,150/255,100/255]},{Thick,RGBColor[0,0,160/255]}},
ImageSize->{chartwidth,400},
AspectRatio->200/chartwidth,
PlotRangePadding->{{3,3},{3,3}},
PlotMarkers->None,
Axes->None,
GridLines->Automatic,
Background->RGBColor[240/255,240/255,243/255],
Frame->True,
RotateLabel-> True,
FrameStyle->{{Directive[Thin,FontFamily->Helvetica,11],Thin},{Directive[Thin,FontFamily->Helvetica,11],Thin}},FrameLabel->{Style["Sweep Number (Decimated Data)",FontSize->11,FontFamily->Helvetica],Style["ADC Value (mV)",FontSize->11,FontFamily->Helvetica],charttitle,""}
];
fore=Rasterize[plot9,ImageSize->chartwidth+110];
mid=ImageCompose[back,fore,{485,150}];
output[[79]]=mid;

(* Create a cross polarization spectrogram                      *)
output[[81]]=Style["Cross Polarization Spectrogram – Decimated    ",Bold];

output[[83]]="         Color file:  "<>xpcolorfile;
output[[84]]="       Color offset:  "<>ToString[xcoloroffset];
output[[85]]="         Color gain:  "<>ToString[PaddedForm[xcolorgain,{2,1},NumberPadding->{"","0"}]];
output[[86]]="               Data:  "<>ToString[sweeps]<>" sweeps";
output[[87]]="       Decimated to:  "<>ToString[chartwidth]<>" sweeps  (decimation factor "<>ToString[PaddedForm[N[shrinkratio],{10,3},NumberPadding->{"","0"}]]<>" )";
output[[88]]="  Decimation method:  "<>"Mean (arithmetic average)";
output[[89]]="            RCP max:  "<>ToString[xchartmax]<>" mV";
output[[90]]="            LCP max:  "<>ToString[xchartmin]<>" mV";

plot5=ArrayPlot[Transpose[dataXadj],
PixelConstrained->False,
ImageSize->{1000,300},
ColorFunctionScaling->False,
ColorFunction->Function[{z},RGBColor[ColorsInTx[[z]]]],
PlotRangePadding->{{10,10},{10,10}},
Background->None,
Frame->True,
FrameTicks->True,
RotateLabel-> True,
FrameStyle->{{Directive[Thin,FontFamily->Helvetica,11],Directive[Thin,FontFamily->Helvetica,11]},{Directive[Thin,FontFamily->Helvetica,11],Directive[Thin,FontFamily->Helvetica,11]}},FrameLabel->{Style["Channel",FontSize->11,FontFamily->Helvetica],Style["",FontSize->11,FontFamily->Helvetica],"Frequency (MHz)",""}
];
back=Image[Raster[Table[{240/255,240/255,243/255},{i,1,275},{j,1,1000}]],ImageSize->{1000,275}];
mid=ImageCompose[back,plot5,{500,138}];
output[[91]]=mid;

(* Create a cross polarization spectrogram                      *)
output[[93]]=Style["Cross Polarization Spectrogram – Decimated – Correction Array Applied    ",Bold];

output[[95]]="         Color file:  "<>xpcolorfile;
output[[96]]="       Color offset:  "<>ToString[xcoloroffset];
output[[97]]="         Color gain:  "<>ToString[PaddedForm[xcolorgain,{2,1},NumberPadding->{"","0"}]];
output[[98]]="               Data:  "<>ToString[sweeps]<>" sweeps";
output[[99]]="       Decimated to:  "<>ToString[chartwidth]<>" sweeps  (decimation factor "<>ToString[PaddedForm[N[shrinkratio],{10,3},NumberPadding->{"","0"}]]<>" )";
output[[100]]="  Decimation method:  "<>"Mean (arithmetic average)";
output[[101]]="            RCP max:  "<>ToString[xchartmaxcorr]<>" mV";
output[[102]]="            LCP max:  "<>ToString[xchartmincorr]<>" mV";

plot10=ArrayPlot[Transpose[dataXadjcorr],
PixelConstrained->False,
ImageSize->{1000,300},
ColorFunctionScaling->False,
ColorFunction->Function[{z},RGBColor[ColorsInTx[[z]]]],
PlotRangePadding->{{10,10},{10,10}},
Background->None,
Frame->True,
FrameTicks->True,
RotateLabel-> True,
FrameStyle->{{Directive[Thin,FontFamily->Helvetica,11],Directive[Thin,FontFamily->Helvetica,11]},{Directive[Thin,FontFamily->Helvetica,11],Directive[Thin,FontFamily->Helvetica,11]}},FrameLabel->{Style["Channel",FontSize->11,FontFamily->Helvetica],Style["",FontSize->11,FontFamily->Helvetica],"Frequency (MHz)",""}
];
back=Image[Raster[Table[{240/255,240/255,243/255},{i,1,275},{j,1,1000}]],ImageSize->{1000,275}];
mid=ImageCompose[back,plot10,{500,138}];
output[[103]]=mid;

(* Output the results                                           *)
g=Panel[
Column[{
Column[output[[10;;11]]],
Column[output[[12;;52]],Alignment->":"],
output[[54]],
Column[output[[55;;56]]],
Column[output[[57;;65]],Alignment->":"],
output[[67]],
Column[output[[68;;70]]],
Column[output[[71;;73]]],
Column[output[[74;;76]]],
Column[output[[77;;79]]],
Column[output[[80;;82]]],
Column[output[[83;;90]],Alignment->":"],
output[[91]],
Column[output[[92;;94]]],
Column[output[[95;;102]],Alignment->":"],
output[[103]]
}],
Column[output[[1;;3]]],Background->RGBColor[240/255,240/255,243/255]
];
g
Export[rptfilename,g]
Export[spgmfilename,output[[54]]]
Export[spgmfilenamecorr,output[[67]]]
Export[spgmAfilename,plot1]
Export[spgmBfilename,plot2]
Export[spgmAfilenamecorr,plot6]
Export[spgmBfilenamecorr,plot7]
Export[contfilename,output[[70]]]
Export[contfilenamecorr,output[[73]]]
Export[stripfilename,output[[76]]]
Export[stripfilenamecorr,output[[79]]]
Export[spgmXfilename,output[[91]]]
Export[spgmXfilenamecorr,output[[103]]]



(* 

Chart Background 240,240,243
Chart Notes, Legend 245,245,248

Add temperature color scale cal'd in K on one side and Jy on the other.

Do data decimation FIRST to keep array size small as possible.

Put data limiting function ahead of colorization routine - may be able to get rid of colorization routine altogether since using color function in arrayplot

ImageCrop or ImageTake to get rid of too-wide borders on dual spectrogram

Instead of making raster arrays containing the color triplets,
Try ArrayPlot[] with ColorFunction-->ColorsInT[[data point value]]
eg. ColorFunction->Function[{x,y,z},Hue[z]]
ColorFunction->Function[{a},ColorsInT[[a]]]
The use of ArrayPlot[] will allow the use of tick marks, etc to do axis labeling 
or
ColorFunction->Function[{a},RGBColor[ColorsInT[[a,1]],ColorsIn[[a,2]],ColorsIn[[a,3]]]]

try to save pdf using lossless compression on images

*)

DPS Data Analysis - OLD  - Full analysis of DPS data

(********************************************************************************************)
(*                                                                                          *)
(*  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.                     *)
(*                                                                                          *)
(********************************************************************************************)

(* Set file name of input file                                  *)
filename="R:\\Observation Records\\AJ4CO_DPS-XTFD\\1305\\130518\\130518164600.sps";

(* Set color offset and color gain for RCP and LCP spectrograms *)
coloroffset=1500;
colorgain=1.2;

(* Set Xpol spectrogram color gain and color offset             *)
xcolorgain=3;
xcoloroffset=-2048;

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

(* Set color file path and name for RCP & LCP spectrograms      *)
cpcolorfile="C:\\Program Files\\Spectrograph\\AJ4CO-Rainbow.txt";

(* Set color file path and name for X-polarization spectrogram  *)
xpcolorfile="C:\\Program Files\\Spectrograph\\CrossPolarization.txt";

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

(* Set strip chart channel number *)
stripchn=110;

(* Create output filenames *)
rptfilename=StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>"_spectrogram_inspection.pdf";
spgmfilename=StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>"_spectrogram_DCP.png";
spgmAfilename=StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>"_spectrogram_RCP.png";
spgmBfilename=StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>"_spectrogram_LCP.png";
spgmXfilename=StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>"_spectrogram_X-pol.png";
contfilename=StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>"_continua.png";
stripfilename=StringTake[filename,StringPosition[filename,".sps"][[1,1]]-1]<>"_strip_chart_channel_"<>ToString[stripchn]<>".png";

(* 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;
(* Graphics[Raster[{ColorsInT}],AspectRatio->.2] *)

(* Repeat above for Xpol color file                             *)
xColorInFile=Import[xpcolorfile,"Data"];
xcolorlo=1;
xcolorhi=Dimensions[xColorInFile][[1]];
xColorsIn ={
FractionalPart[FractionalPart[xColorInFile[[All,1]]/65536]*65536/256]*256,
IntegerPart[FractionalPart[xColorInFile[[All,1]]/65536]*65536/256],
IntegerPart[xColorInFile[[All,1]]/65536]
};
xColorsInT=Transpose[xColorsIn];
ColorsInTx=xColorsInT/255;
(* Graphics[Raster[{ColorsInTx}],AspectRatio->.2] *)

(* Open one SPS file for input                                  *)
infile=OpenRead[filename,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 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=Round[(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=Round[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 >"];

(* Create temp table and read in spectro data.                  *)
(* 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)                *)
dataAdecimated=Table[0,{i,1,chartwidth},{j,1,chartchannels}];
dataBdecimated=Table[0,{i,1,chartwidth},{j,1,chartchannels}];
(* Create input record accumulators to sum & average input data *)
indataArec=Table[0,{i,1,chartchannels}];
indataBrec=Table[0,{i,1,chartchannels}];
shrinkratio=sweeps/chartwidth;
(* Initialize raw data min-max accumulators                     *)
chartAmaxR=0;
chartAminR=2^nbits-1;
chartBmaxR=0;
chartBminR=2^nbits-1;
(* Read data for all ouput chart records                        *)
Do[
{
(* Write averaged input records to one chart record          *)
Do[
{
(* Read input data record                                  *)
Do[
{
(* Read Side B data (side B comes first in the SPS data) *)
tempB=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
chartBmaxR=Max[chartBmaxR,tempB];
chartBminR=Min[chartBminR,tempB];
indataBrec[[i]]=indataBrec[[i]]+tempB;
(* Read Side A data                                      *)
tempA=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
chartAmaxR=Max[chartAmaxR,tempA];
chartAminR=Min[chartAminR,tempA];
indataArec[[i]]=indataArec[[i]]+tempA;
}
(* cycle over 2 x number of channels *)
,{i,1,chartchannels}
];
(* Read EOR marker                                         *)
dummy=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
}
(* cycle over number of input sweeps to be averaged         *)
,{j,Round[1+(z-1)*(shrinkratio)],Round[z*shrinkratio]} 
];
(*  RSS and DPS data have channel 200 rolled over above channel 1. *)
(*  These two lines put channel 200 data back where it belongs.    *)
indataArec=RotateLeft[indataArec];
indataBrec=RotateLeft[indataBrec];
(* Average summed input records and write to decimated array *)
(* See "Subsampling Area Averaging Method.xls" for details   *)
dataAdecimated[[z]]=Round[(indataArec/(Round[z*shrinkratio]-Round[1+(z-1)*shrinkratio]+1))];
dataBdecimated[[z]]=Round[(indataBrec/(Round[z*shrinkratio]-Round[1+(z-1)*shrinkratio]+1))];
(* Zero out input record accumulators before next round of averaging *)
indataArec=indataArec*0;
indataBrec=indataBrec*0;
}
(* cycle over output chart width in pixels, sweeps, records, whatever *)
,{z,1,chartwidth}
];
Close[infile];

(* Find min and max of decimated data over all data points      *)
chartAmaxD=Round[Max[dataAdecimated]];
chartAminD=Round[Min[dataAdecimated]];
chartBmaxD=Round[Max[dataBdecimated]];
chartBminD=Round[Min[dataBdecimated]];

(* Apply color offset and color gain to raw data, just like RSS *) 
dataAadj=Round[(dataAdecimated-coloroffset)*colorgain];
dataBadj=Round[(dataBdecimated-coloroffset)*colorgain];

(* Coerce adjusted data that lays outside n-bit range          *)
datalo=1;
datahi=2^nbits-1;  
(* Coerce Side A data                                          *)
Do[{
If[
dataAadj[[i,j]]≤datahi,
{If[
dataAadj[[i,j]]≥datalo,
{dataAadj[[i,j]]=dataAadj[[i,j]]},
{dataAadj[[i,j]]=datalo}
]
},
{dataAadj[[i,j]]=datahi}
]
}
,{i,1,chartwidth}
,{j,1,chartchannels}
]
(* Coerce Side B data                                          *)
Do[{
If[
dataBadj[[i,j]]≤datahi,
{If[
dataBadj[[i,j]]≥datalo,
{dataBadj[[i,j]]=dataBadj[[i,j]]},
{dataBadj[[i,j]]=datalo}
]
},
{dataBadj[[i,j]]=datahi}
]
}
,{i,1,chartwidth}
,{j,1,chartchannels}
]

(* Create X-pol data array                                     *)
dataXdecimated=dataAdecimated;
dataXdecimated=dataAdecimated-dataBdecimated;
xchartmax=Round[Max[dataXdecimated]];
xchartmin=Round[Min[dataXdecimated]];
dataXadj=Round[dataXdecimated*xcolorgain-xcoloroffset];
(* Coerce adjusted data that lays outside n-bit range          *)
datalo=1;
datahi=2^nbits-1;  
(* Coerce X-pol data                                           *)
Do[{
If[
dataXadj[[i,j]]≤datahi,
{If[
dataXadj[[i,j]]≥datalo,
{dataXadj[[i,j]]=dataXadj[[i,j]]},
{dataXadj[[i,j]]=datalo}
]
},
{dataXadj[[i,j]]=datahi}
]
}
,{i,1,chartwidth}
,{j,1,chartchannels}
]

(* Calculate time averaged continua for Side A & Side B data    *)
dataAM=N[Mean[dataAdecimated]];
dataBM=N[Mean[dataBdecimated]];

(* Define galactic background model, output in kK               *)
GB[f_]:=1.07*10^5*f^-2.56;
(* Define function to convert GB kK to ADC mV                   *)
GBmV[f_]:=10 Log[10,GB[f]/GB[32.92]]*25+1670;
(* Create table of GB model ADC values in mV                    *)
dataGB=Table[GBmV[f],{f,33-0.08,17,-0.08}];

(* Print formatted file header data to screen                   *)
output=Table["",{i,1,150}];
output[[1]]=Style["DPS Data Analysis",Larger,Darker[Darker[Blue]],Bold];
output[[2]]="Run date: "<>DateString[{"Day","-","MonthNameShort","-","Year"," ","Hour24",":","Minute"," Local"}];

output[[10]]=Style["File Information – Raw Data",Bold];

output[[12]]="         Local name:  "<>chartlocalname;
output[[13]]="         Input file:  "<> filename;
output[[14]]="        RSS version:  "<>chartversionF;
output[[15]]="        Chart start:  "<>chartstartF;
output[[16]]="          Chart end:  "<>chartendF;

output[[18]]="         Local name:  "<>chartlocalname;
output[[19]]="     Chart location:  "<>chartlocation;
output[[20]]="           Latitude:  "<>ToString[PaddedForm[chartlat,{8,6},NumberPadding->{"0","0"},NumberSigns->{"(S) -","(N) +"}]]<>" degrees";
output[[21]]="          Longitude:  "<>ToString[PaddedForm[chartlong,{8,6},NumberPadding->{"0","0"},NumberSigns->{"(W) -","(E) +"}]]<>" degrees";
output[[22]]="          Time zone:  UTC+"<>ToString[chartTZ];
output[[23]]="           Observer:  "<>chartauthor;

output[[25]]="        Data source:  "<>chartsource;
output[[26]]="      Data duration:  "<>ToString[PaddedForm[chartduration,{10,3},NumberPadding->{"","0"}]]<>" s";
output[[27]]="        Sweep start:  "<>ToString[PaddedForm[freqhi,{10,3},NumberPadding->{"","0"}]]<>" MHz";
output[[28]]="          Sweep end:  "<>ToString[PaddedForm[freqlo,{10,3},NumberPadding->{"","0"}]]<>" MHz";
output[[29]]="             Sweeps:  "<>ToString[sweeps];
output[[30]]="         Sweep rate:  "<>ToString[PaddedForm[chartsweeprate,{6,3},NumberPadding->{"","0"}]]<>" Hz";
output[[31]]="           Channels:  "<>ToString[chartchannels];
output[[32]]=" Channel separation:  "<>ToString[PaddedForm[channelsep,{6,3},NumberPadding->{"","0"}]]<>" kHz";
output[[33]]=" Channel dwell time:  "<>ToString[channeldwelltime]<>" μs";

output[[35]]="RCP ADC range (raw):  "<>ToString[chartAminR]<>" mV min  to  "<>ToString[chartAmaxR]<>" mV max";
output[[36]]="LCP ADC range (raw):  "<>ToString[chartBminR]<>" mV min  to  "<>ToString[chartBmaxR]<>" mV max";

output[[38]]="   Data note length:  "<>ToString[chartnotelength];
output[[39]]="          Data note:  "<>chartnote;

output[[41]]=Style["Spectrograms – Decimated    ",Bold];

output[[43]]="         Color file:  "<>cpcolorfile;
output[[44]]="       Color offset:  "<>ToString[coloroffset];
output[[45]]="         Color gain:  "<>ToString[PaddedForm[colorgain,{2,1},NumberPadding->{"","0"}]];
output[[46]]="               Data:  "<>ToString[sweeps]<>" sweeps";
output[[47]]="       Decimated to:  "<>ToString[chartwidth]<>" sweeps  (decimation factor "<>ToString[PaddedForm[N[shrinkratio],{10,3},NumberPadding->{"","0"}]]<>" )";
output[[48]]="  Decimation method:  "<>"Mean (arithmetic average)";
output[[49]]="RCP ADC range (decimated):  "<>ToString[chartAminD]<>" mV min  to  "<>ToString[chartAmaxD]<>" mV max";
output[[50]]="LCP ADC range (decimated):  "<>ToString[chartBminD]<>" mV min  to  "<>ToString[chartBmaxD]<>" mV max";

plot1=ArrayPlot[Transpose[dataAadj],
PixelConstrained->False,
ImageSize->{1000,300},
ColorFunctionScaling->False,
ColorFunction->Function[{z},RGBColor[ColorsInT[[z]]]],
PlotRangePadding->{{10,10},{10,10}},
Background->None,
Frame->True,
FrameTicks->True,
RotateLabel-> True,
FrameStyle->{{Directive[Thin,FontFamily->Helvetica,11],Directive[Thin,FontFamily->Helvetica,11]},{Directive[Thin,FontFamily->Helvetica,11],Directive[Thin,FontFamily->Helvetica,11]}},FrameLabel->{Style["RCP   Channel",FontSize->11,FontFamily->Helvetica],Style["",FontSize->11,FontFamily->Helvetica],"RCP   Frequency (MHz)",""}
];

plot2=ArrayPlot[Transpose[dataBadj],
PixelConstrained->False,
ImageSize->{1000,300},
ColorFunctionScaling->False,
ColorFunction->Function[{z},RGBColor[ColorsInT[[z]]]],
PlotRangePadding->{{10,10},{10,10}},
Background->None,
Frame->True,
FrameTicks->True,
RotateLabel-> True,
FrameStyle->{{Directive[Thin,FontFamily->Helvetica,11],Directive[Thin,FontFamily->Helvetica,11]},{Directive[Thin,FontFamily->Helvetica,11],Directive[Thin,FontFamily->Helvetica,11]}},FrameLabel->{Style["LCP   Channel",FontSize->11,FontFamily->Helvetica],Style["",FontSize->11,FontFamily->Helvetica],"LCP   Frequency (MHz)",""}
];
back=Image[Raster[Table[{240/255,240/255,243/255},{i,1,546},{j,1,1000}]],ImageSize->{1000,546}];
mid=ImageCompose[back,plot1,{500,396}];
dual=ImageCompose[mid,plot2,{500,150}];
output[[75]]=dual;
output[[80]]=Style["Continua – Decimated",Bold];

(* Create a continuum plot                                      *)
charttitle=Column[{"",
Style["RCP & LCP Continua – average of "<>ToString[chartwidth]<>" sweeps"<>" decimated from "<>ToString[sweeps]<>" sweeps  ("<>ToString[Round[chartduration]]<>" seconds in original data)",Plain,FontSize->12,FontFamily->Helvetica],
Row[{
Style["RCP = Green",Plain,FontSize->12,FontFamily->Helvetica,RGBColor[0,150/255,100/255]],
Style["LCP = RED",Plain,FontSize->12,FontFamily->Helvetica,RGBColor[190/255,39/255,39/255]],
Style["GB Model = Blue (for slope ONLY)",Plain,FontSize->12,FontFamily->Helvetica,RGBColor[0,0,160/255]]
},"        "]
},Center];
plot3=ListLinePlot[{dataBM,dataAM,dataGB},
DataRange->{1,chartchannels},
PlotRange->{1500,2100},
PlotStyle->{{Thick,RGBColor[190/255,39/255,39/255]},{Thick,RGBColor[0,150/255,100/255]},{Thick,RGBColor[0,0,160/255]}},
AspectRatio->200/chartwidth,
PlotRangePadding->{{3,3},{3,3}},
PlotMarkers->None,
Axes->None,
GridLines->Automatic,
Background->RGBColor[240/255,240/255,243/255],
Frame->True,
RotateLabel-> True,
FrameStyle->{{Directive[Thin,FontFamily->Helvetica,11],Thin},{Directive[Thin,FontFamily->Helvetica,11],Thin}},FrameLabel->{Style["Channel Number",FontSize->11,FontFamily->Helvetica],Style["Average ADC Value (mV)",FontSize->11,FontFamily->Helvetica],charttitle,""}
];
back=Image[Raster[Table[{240/255,240/255,243/255},{i,1,300},{j,1,1000}]],ImageSize->{1000,300}];
fore=Rasterize[plot3,ImageSize->chartwidth+110];
mid=ImageCompose[back,fore,{485,150}];
output[[81]]=mid;

output[[85]]=Style["Strip Chart – Decimated",Bold];
(* Create a continuum plot                                      *)
charttitle=Column[{"",
Style["RCP & LCP Strip Chart – "<>ToString[chartwidth]<>" sweeps"<>" decimated from "<>ToString[sweeps]<>" sweeps – Channel "<>ToString[stripchn],Plain,FontSize->12,FontFamily->Helvetica],
Row[{
Style["RCP = Green",Plain,FontSize->12,FontFamily->Helvetica,RGBColor[0,150/255,100/255]],
Style["LCP = RED",Plain,FontSize->12,FontFamily->Helvetica,RGBColor[190/255,39/255,39/255]],
Style["GB Model = Blue (NOT YET MEANINGFUL)",Plain,FontSize->12,FontFamily->Helvetica,RGBColor[0,0,160/255]]
},"        "]
},Center];
plot4=ListLinePlot[{dataBdecimated[[All,stripchn]],dataAdecimated[[All,stripchn]],Table[dataGB[[stripchn]],{i,1,chartwidth}]},
DataRange->{1,chartwidth},
PlotRange->{1500,2500},
PlotStyle->{{Thick,RGBColor[190/255,39/255,39/255]},{Thick,RGBColor[0,150/255,100/255]},{Thick,RGBColor[0,0,160/255]}},
ImageSize->{chartwidth,400},
AspectRatio->200/chartwidth,
PlotRangePadding->{{3,3},{3,3}},
PlotMarkers->None,
Axes->None,
GridLines->Automatic,
Background->RGBColor[240/255,240/255,243/255],
Frame->True,
RotateLabel-> True,
FrameStyle->{{Directive[Thin,FontFamily->Helvetica,11],Thin},{Directive[Thin,FontFamily->Helvetica,11],Thin}},FrameLabel->{Style["Sweep Number (Decimated Data)",FontSize->11,FontFamily->Helvetica],Style["ADC Value (mV)",FontSize->11,FontFamily->Helvetica],charttitle,""}
];
fore=Rasterize[plot4,ImageSize->chartwidth+110];
mid=ImageCompose[back,fore,{485,150}];
output[[86]]=mid;

(* Create a cross polarization spectrogram                      *)
output[[90]]=Style["Cross Polarization Spectrogram – Decimated    ",Bold];

output[[92]]="         Color file:  "<>xpcolorfile;
output[[93]]="       Color offset:  "<>ToString[xcoloroffset];
output[[94]]="         Color gain:  "<>ToString[PaddedForm[xcolorgain,{2,1},NumberPadding->{"","0"}]];
output[[95]]="               Data:  "<>ToString[sweeps]<>" sweeps";
output[[96]]="       Decimated to:  "<>ToString[chartwidth]<>" sweeps  (decimation factor "<>ToString[PaddedForm[N[shrinkratio],{10,3},NumberPadding->{"","0"}]]<>" )";
output[[97]]="  Decimation method:  "<>"Mean (arithmetic average)";
output[[98]]="            RCP max:  "<>ToString[xchartmax]<>" mV";
output[[99]]="            LCP max:  "<>ToString[xchartmin]<>" mV";

plot5=ArrayPlot[Transpose[dataXadj],
PixelConstrained->False,
ImageSize->{1000,300},
ColorFunctionScaling->False,
ColorFunction->Function[{z},RGBColor[ColorsInTx[[z]]]],
PlotRangePadding->{{10,10},{10,10}},
Background->None,
Frame->True,
FrameTicks->True,
RotateLabel-> True,
FrameStyle->{{Directive[Thin,FontFamily->Helvetica,11],Directive[Thin,FontFamily->Helvetica,11]},{Directive[Thin,FontFamily->Helvetica,11],Directive[Thin,FontFamily->Helvetica,11]}},FrameLabel->{Style["Channel",FontSize->11,FontFamily->Helvetica],Style["",FontSize->11,FontFamily->Helvetica],"Frequency (MHz)",""}
];
back=Image[Raster[Table[{240/255,240/255,243/255},{i,1,275},{j,1,1000}]],ImageSize->{1000,275}];
mid=ImageCompose[back,plot5,{500,150}];
output[[100]]=mid;

(* Output the results                                           *)
g=Panel[
Column[{
Column[output[[10;;11]]],
Column[output[[12;;50]],Alignment->":"],
output[[75]],
Column[output[[80;;81]]],
Column[output[[84;;86]]],
Column[output[[89;;90]]],
Column[output[[91;;99]],Alignment->":"],
output[[100]]
}],
Column[output[[1;;3]]],Background->RGBColor[240/255,240/255,243/255]
];
g
Export[rptfilename,g]
Export[spgmfilename,dual]
Export[spgmAfilename,plot1]
Export[spgmBfilename,plot2]
Export[contfilename,output[[81]]]
Export[stripfilename,output[[86]]]
Export[spgmXfilename,output[[100]]]









(* 

Chart Background 240,240,243
Chart Notes, Legend 245,245,248

Add temperature color scale cal'd in K on one side and Jy on the other.

Do data decimation FIRST to keep array size small as possible.

Put data limiting function ahead of colorization routine - may be able to get rid of colorization routine altogether since using color function in arrayplot

ImageCrop or ImageTake to get rid of too-wide borders on dual spectrogram

Instead of making raster arrays containing the color triplets,
Try ArrayPlot[] with ColorFunction-->ColorsInT[[data point value]]
eg. ColorFunction->Function[{x,y,z},Hue[z]]
ColorFunction->Function[{a},ColorsInT[[a]]]
The use of ArrayPlot[] will allow the use of tick marks, etc to do axis labeling 
or
ColorFunction->Function[{a},RGBColor[ColorsInT[[a,1]],ColorsIn[[a,2]],ColorsIn[[a,3]]]]

try to save pdf using lossless compression on images

*)

(********************************************************************************************)
(*                                                                                          *)
(*  Load and apply correction files to raw spectro data, display corrected spectrograms     *)
(*                                                                                          *)
(*                                                                                          *)
(********************************************************************************************)

(* Load DPS Side A correction array & transpose since arrays    *)
(* exported and subsequently imported undergo transposition.    *) 
dataACorr=dataA;
corrarrayA=Transpose[Import[corrarrayAinfile]][[1]];
(* Apply correction array to Side A data                        *)
Do[dataACorr[[i]]=dataA[[i,All]]+corrarrayA,{i,1,chartrecs}];
(* Load Side B correction array & apply to Side B data          *) 
dataBCorr=dataB;
corrarrayB=Transpose[Import[corrarrayBinfile]][[1]];
Do[dataBCorr[[i]]=dataB[[i,All]]+corrarrayB,{i,1,chartrecs}];

(* Apply color offset and color gain to corrected data          *) 
dataAadj=Round[(dataACorr-coloroffset)*colorgain];
dataBadj=Round[(dataBCorr-coloroffset)*colorgain];

Do[{
If[
dataAadj[[i,j]]<=4096,
{If[
dataAadj[[i,j]]≥1,
{dataAadj[[i,j]]=ColorsInT[[dataAadj[[i,j]]]]},
{dataAadj[[i,j]]=ColorsInT[[1]]}
]
},
{dataAadj[[i,j]]=ColorsInT[[4096]]}
]
}
,{i,1,chartrecs}
,{j,1,chartchannels}
]

Do[{
If[
dataBadj[[i,j]]<=4096,
{If[
dataBadj[[i,j]]≥1,
{dataBadj[[i,j]]=ColorsInT[[dataBadj[[i,j]]]]},
{dataBadj[[i,j]]=ColorsInT[[1]]}
]
},
{dataBadj[[i,j]]=ColorsInT[[4096]]}
]
}
,{i,1,chartrecs}
,{j,1,chartchannels}
]

DPS File Reviewer - With Cross Polarization Analysis - OLD

(* Open one file for input *)
filename="R:\\Observation Records\\AJ4CO_DPS-XTFD\\1305\\130522\\130522130700.sps";

(* Set input file path and name for DPS Xpol SPS file           *)
outfilename="R:\\Observation Records\\AJ4CO_DPS-XTFD\\1305\\130510\\130510181100_Xpol.sps";

(* Set color file path and name for RCP & LCP spectrograms      *)
cpcolorfile="C:\\Program Files\\Spectrograph\\AJ4CO-Rainbow.txt";

(* Set color file path and name for Xpol spectrogram            *)
xpolcolorfile="C:\\Program Files\\Spectrograph\\CrossPolarization.txt";

(* Set correction file paths and names                          *)
corrarrayAinfile="R:\\Observation Records\\AJ4CO_DPS-XTFD\\1305\\130510\\corrarrayA.csv";
corrarrayBinfile="R:\\Observation Records\\AJ4CO_DPS-XTFD\\1305\\130510\\corrarrayB.csv";

(* Set output PNG file path and name                            *)
outpngfile="R:\\Observation Records\\AJ4CO_DPS-XTFD\\1305\\130510\\130510181100_Xpol.png";

(* Set color offset and color gain for RCP and LCP spectrograms *)
coloroffset=1500;
colorgain=1.2;

(* Set Xpol spectrogram color gain and color offset             *)
xpolcgain=3;
xpolcoffset=2048;

(* Set output spectrogram widths in pixels *)
spectrowidth=900;


(********************************************************************************************)
(*                                                                                          *)
(*  Load file, read header and raw data, display data parameters and raw spectrographs      *)
(*                                                                                          *)
(*                                                                                          *)
(********************************************************************************************)

infile=OpenRead[filename,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 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=Round[(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=(freqhi-freqlo)/chartchannels*10^3;
channeldwelltime=Round[chartduration/(chartrecs*chartchannels)*10^6];
sweeps=ToExpression[StringTake[chartnote,{ StringPosition[chartnote,"SWEEPS"][[1,2]]+1 ,StringPosition[chartnote,"ÿLOWF"][[1,1]]-1}]];

(* Create temp table and read in spectro data.                  *)
(* 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)                *)
indata=Table[0,{i,1,sweeps},{j,1,chartchannels*2+1}];
Do[{
indata[[i,j]]=BinaryRead[infile,"UnsignedInteger16",ByteOrdering->+1];
}
,{i,1,sweeps}
,{j,1,chartchannels*2+1}
];
Close[infile];

(* Verify number of records present in data file by checking    *)
(* count of EOR markers.                                        *)
(* The EOR marker is two consecutive 0xFE characters, ASCII "þ" *)
(* or 65278 decomal when read as a 2-byte unsigned integer.     *)
chartrecs=Count[indata,65278,2];

chartrecs
sweeps

(*
(* Calculate chart sweep rate                                   *)
chartsweeprate=chartrecs/chartduration;

(* Find min and max of raw data over all data points            *)
chartmaxF=Max[indata[[1;;chartrecs,1;;chartchannels*2]]];
chartminF=Min[indata[[1;;chartrecs,1;;chartchannels*2]]];

(* Print formatted file header data to screen                   *)
Print["DPS Cross Polarization Analysis"]
Print["Analysis run date: "<>DateString[{"Day","-","MonthNameShort","-","Year"," ","Hour24",":","Minute"," Local"}]]
Print[]
Print["         Local name: ",chartlocalname]
Print["         Input file: ", filename]
Print["        RSS version: ",chartversionF]
Print["        Chart start: ",chartstartF]
Print["          Chart end: ",chartendF]
Print[]
Print["         Local name: ",chartlocalname]
Print["     Chart location: ",chartlocation]
Print["           Latitude: ",ToString[PaddedForm[chartlat,{8,6},NumberPadding->{"0","0"},NumberSigns->{"(S) -","(N) +"}]]<>" degrees"]
Print["          Longitude: ",ToString[PaddedForm[chartlong,{8,6},NumberPadding->{"0","0"},NumberSigns->{"(W) -","(E) +"}]]<>" degrees"]
Print["          Time zone: UTC+",chartTZ]
Print["           Observer: ",chartauthor]
Print[]
Print["        Data Source: ",chartsource]
Print["      Data duration: ",ToString[PaddedForm[chartduration,{10,3},NumberPadding->{"","0"}]]," s"]
Print["        Sweep Start: ",ToString[PaddedForm[freqhi,{10,6},NumberPadding->{"","0"}]]," MHz"]
Print["          Sweep End: ",ToString[PaddedForm[freqlo,{10,6},NumberPadding->{"","0"}]]," MHz"]
Print["             Sweeps: ",chartrecs]
Print["         Sweep rate: ",ToString[PaddedForm[chartsweeprate,{6,3},NumberPadding->{"","0"}]]," Hz"]
Print["           Channels: ",chartchannels]
Print[" Channel separation: ",ToString[PaddedForm[channelsep,{6,3},NumberPadding->{"","0"}]]," kHz"]
Print[" Channel dwell time: ",ToString[channeldwelltime]<>" μs"]
Print[]
Print["           Data Max: ",ToString[chartmaxF]<>" mV"]
Print["           Data Min: ",ToString[chartminF]<>" mV"]
Print[]
Print["  Chart note length: ",chartnotelength]
Print["       Chart note 1: ",chartnote]

(* Create 12-bit data from input file array; this strips out    *)
(* the EOR markers from the array.                              *)
(* DPS SPS files are already 12-bit, so no multiplication is    *)
(* necessary.                                                   *)
indata12bit=indata[[1;;chartrecs,1;;chartchannels*2]];

(* Create DPS Side A nd Side B arrays.                          *)
(* Normal operation: Side A = RCP, Side B = LCP                 *)
(* Extract RCP and LCP data from interleaved raw data           *)
(* Raw data record format:                                      *)
(* Side B1, Side A1, Side B2, Side A2 ... Side Bn, Side An, EOR *)
dataA=Table[0.0,{i,1,chartrecs},{j,1,chartchannels}];
dataB=Table[0.0,{i,1,chartrecs},{j,1,chartchannels}];
Do[
Do[
{dataA[[i,j]]=indata12bit[[i,j*2]],
dataB[[i,j]]=indata12bit[[i,j*2-1]]}
,{j,1,chartchannels}
],{i,1,chartrecs}
];

(* Apply color offset and color gain to raw data, just like RSS *) 
dataAadj=Round[(dataA-coloroffset)*colorgain];
dataBadj=Round[(dataB-coloroffset)*colorgain];

(* ***************************************************************************************************************)
(* insert data shrinkage routine here *)
(* this is the smoothed data that will be the source data for all subsequent operations *)
dataAshrunk=dataAadj;
dataBshrunk=dataBadj;
chartshrunkrecs=chartrecs;

(* Open and read color file data for circular pol spectrograms  *)
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;

(*
Graphics[Raster[{ColorsInT}],AspectRatio->.1]
*)

ts=ImageRotate[Image[Graphics[Raster[{ColorsInT}],AspectRatio->.1]],π/2]


(* Create raster arrays of RGB triplets of raw spectro data     *)
dataAshrunkraster=dataAshrunk;
dataBshrunkraster=dataBshrunk;
(* Replace raw data values with color triplets based on color   *)
(* file data                                                    *)
(* Colorize DPS Side A shrunk data                              *)  
Do[{
If[
dataAshrunkraster[[i,j]]≤colorhi,
{If[
dataAshrunkraster[[i,j]]≥colorlo,
{dataAshrunkraster[[i,j]]=ColorsInT[[dataAshrunkraster[[i,j]]]]},
{dataAshrunkraster[[i,j]]=ColorsInT[[colorlo]]}
]
},
{dataAshrunkraster[[i,j]]=ColorsInT[[colorhi]]}
]
}
,{i,1,chartshrunkrecs}
,{j,1,chartchannels}
]
(* Colorize DPS Side B shrunk data                              *)
Do[{
If[
dataBshrunkraster[[i,j]]≤colorhi,
{If[
dataBshrunkraster[[i,j]]≥colorlo,
{dataBshrunkraster[[i,j]]=ColorsInT[[dataBshrunkraster[[i,j]]]]},
{dataBshrunkraster[[i,j]]=ColorsInT[[colorlo]]}
]
},
{dataBshrunkraster[[i,j]]=ColorsInT[[colorhi]]}
]
}
,{i,1,chartshrunkrecs}
,{j,1,chartchannels}
]
(* Transpose arrays to match Mathematica's default for raster   *)
(* images.                                                      *)
dataAshrunkrasterT=Reverse[Transpose[dataAshrunkraster],1];
dataBshrunkrasterT=Reverse[Transpose[dataBshrunkraster],1];

(* Print text info and spectrograms                             *)
Print[]
Print[]
Print["Raw Data Analysis"]
Print[]
Print["         Color File: "<>cpcolorfile]
Print["       Color Offset: ",IntegerString[coloroffset]]
Print["         Color Gain: ",ToString[PaddedForm[colorgain,{2,1},NumberPadding->{"","0"}]]]
Print["       Decimated to: ",IntegerString[chartshrunkrecs]," sweeps"]
Print[]
Print["       Polarization: RCP - raw data"]
ImageResize[Raster[dataAshrunkrasterT, Automatic],{900,200},Resampling-> "Bilinear"]
Print["       Polarization: LCP - raw data"]
ImageResize[Raster[dataBshrunkrasterT, Automatic],{900,200},Resampling-> "Bilinear"]


(********************************************************************************************)
(*                                                                                          *)
(*  Load and apply correction files to raw spectro data, display corrected spectrograms     *)
(*                                                                                          *)
(*                                                                                          *)
(********************************************************************************************)

(* Load DPS Side A correction array & transpose since arrays    *)
(* exported and subsequently imported undergo transposition.    *) 
dataACorr=dataA;
corrarrayA=Transpose[Import[corrarrayAinfile]][[1]];
(* Apply correction array to Side A data                        *)
Do[dataACorr[[i]]=dataA[[i,All]]+corrarrayA,{i,1,chartrecs}];
(* Load Side B correction array & apply to Side B data          *) 
dataBCorr=dataB;
corrarrayB=Transpose[Import[corrarrayBinfile]][[1]];
Do[dataBCorr[[i]]=dataB[[i,All]]+corrarrayB,{i,1,chartrecs}];

(* Apply color offset and color gain to corrected data          *) 
dataAadj=Round[(dataACorr-coloroffset)*colorgain];
dataBadj=Round[(dataBCorr-coloroffset)*colorgain];

(* ***************************************************************************************************************)
(* insert data shrinkage routine here *)
(* this is the smoothed data that will be the source data for all subsequent operations *)
dataAshrunk=dataAadj;
dataBshrunk=dataBadj;
chartshrunkrecs=chartrecs;

(* Create raster arrays of RGB triplets of corrected data       *)
dataAshrunkraster=dataAshrunk;
dataBshrunkraster=dataBshrunk;
(* Replace corrected data values with color triplets based on   *)
(* color file data.                                             *)
(* Colorize DPS Side A corrected & shrunk data                  *)  
Do[{
If[
dataAshrunkraster[[i,j]]≤colorhi,
{If[
dataAshrunkraster[[i,j]]≥colorlo,
{dataAshrunkraster[[i,j]]=ColorsInT[[dataAshrunkraster[[i,j]]]]},
{dataAshrunkraster[[i,j]]=ColorsInT[[colorlo]]}
]
},
{dataAshrunkraster[[i,j]]=ColorsInT[[colorhi]]}
]
}
,{i,1,chartshrunkrecs}
,{j,1,chartchannels}
]
(* Colorize DPS Side B corrected & shrunk data                  *)
Do[{
If[
dataBshrunkraster[[i,j]]≤colorhi,
{If[
dataBshrunkraster[[i,j]]≥colorlo,
{dataBshrunkraster[[i,j]]=ColorsInT[[dataBshrunkraster[[i,j]]]]},
{dataBshrunkraster[[i,j]]=ColorsInT[[colorlo]]}
]
},
{dataBshrunkraster[[i,j]]=ColorsInT[[colorhi]]}
]
}
,{i,1,chartshrunkrecs}
,{j,1,chartchannels}
]
(* Transpose arrays to match Mathematica's default for raster   *)
(* images.                                                      *)
dataAshrunkrasterT=Reverse[Transpose[dataAshrunkraster],1];
dataBshrunkrasterT=Reverse[Transpose[dataBshrunkraster],1];

(* Print text info and spectrograms                             *)
Print[]
Print[]
Print["Data Analysis, Correction Arrays Applied"]
Print[]
Print["         Color File: "<>cpcolorfile]
Print["       Color Offset: ",IntegerString[coloroffset]]
Print["         Color Gain: ",ToString[PaddedForm[colorgain,{2,1},NumberPadding->{"","0"}]]]
Print["       Decimated to: ",IntegerString[chartshrunkrecs]," sweeps"]
Print[]
Print["       Polarization: RCP - raw data"]
ImageResize[Raster[dataAshrunkrasterT, Automatic],{900,200},Resampling-> "Bilinear"]
Print["       Polarization: LCP - raw data"]
ImageResize[Raster[dataBshrunkrasterT, Automatic],{900,200},Resampling-> "Bilinear"]


(********************************************************************************************)
(*                                                                                          *)
(*  Create cross polarization data, assign color triplets, display image, and export        *)
(*  a DPS data file containing Xpol data.                                                   *)
(*                                                                                          *)
(********************************************************************************************)

(* Open and read color file data for Xpol spectrograms  *)
ColorInFile=Import[xpolcolorfile,"Data"];
(* first element of color data                                  *)
lo=1;
(* number of color data levels, usually 4096                    *)
hi=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;

(* Create Xpol data from corrected raw data                     *)
dataCadj=Round[(dataACorr-dataBCorr)*xpolcgain]+xpolcoffset;

(* ***************************************************************************************************************)
(* insert data shrinkage routine here *)
(* this is the smoothed data that will be the source data for all subsequent operations *)
dataCshrunk=dataCadj;
chartshrunkrecs=chartrecs;

(* Create raster arrays of RGB triplets of corrected data       *)
dataCshrunkraster=dataCshrunk;

(* Replace corrected data values with color triplets based on   *)
(* color file data.                                             *)
(* Colorize corrected, X-polarized, & shrunk data               *)  
Do[{
If[
dataCshrunkraster[[i,j]]≤colorhi,
{If[
dataCshrunkraster[[i,j]]≥colorlo,
{dataCshrunkraster[[i,j]]=ColorsInT[[dataCshrunkraster[[i,j]]]]},
{dataCshrunkraster[[i,j]]=ColorsInT[[colorlo]]}
]
},
{dataCshrunkraster[[i,j]]=ColorsInT[[colorhi]]}
]
}
,{i,1,chartshrunkrecs}
,{j,1,chartchannels}
]
(* Transpose arrays to match Mathematica's default for raster   *)
(* images.                                                      *)
dataCshrunkrasterT=Reverse[Transpose[dataCshrunkraster],1];

(* Print text info and Xpol spectrogram                         *)
Print[]
Print[]
Print["Cross Polarization Analysis"]
Print[]
Print["         Color File: "<>xpolcolorfile]
Print["    Xpol Color Gain: ",IntegerString[xpolcgain]]
Print["  Xpol Color Offset: ",IntegerString[xpolcoffset]]
Print["       Decimated to: ",IntegerString[chartshrunkrecs]," sweeps"]
Print[]
Print["       Polarization: RCP (green) minus LCP (red), after RCP & LCP correction arrays applied to raw data"]
cpol=ImageResize[Raster[dataCshrunkrasterT, Automatic],{900,200},Resampling-> "Bilinear"]

(* Save Xpol spectrogram to PNG file                            *)
Export[outpngfile,cpol];
Print[]
Print["Xpol spectrogram saved as "<>outpngfile]

(* Open new SPS file for output                                 *)
outfile=OpenWrite[outfilename,BinaryFormat->True];

(* Write file header info                                       *)
WriteString[outfile,chartversion];
BinaryWrite[outfile,chartstart,"Real64"];
BinaryWrite[outfile,chartend,"Real64"];
BinaryWrite[outfile,chartlat,"Real64"];
BinaryWrite[outfile,chartlong,"Real64"];
BinaryWrite[outfile,chartmax,"Real64"];
BinaryWrite[outfile,chartmin,"Real64"];
BinaryWrite[outfile,chartTZ,"Integer16"];
WriteString[outfile,chartsource];
WriteString[outfile,chartauthor];
WriteString[outfile,chartlocalname];
WriteString[outfile,chartlocation];
BinaryWrite[outfile,chartchannels,"Integer16"];
BinaryWrite[outfile,chartnotelength,"Integer32"];
WriteString[outfile,chartnote];

(* Make sure data to be written lies in a 12-bit range          *)
Do[{
If[
dataCadj[[i,j]]≤colorhi,
{If[
dataCadj[[i,j]]≥colorlo,
{dataCadj[[i,j]]=dataCadj[[i,j]]},
{dataCadj[[i,j]]=colorlo}
]
},
{dataCadj[[i,j]]=colorhi}
]
}
,{i,1,chartrecs}
,{j,1,chartchannels}
]

(* Write Xpol data using DPS file format                        *)
Do[
{
Do[
{
BinaryWrite[outfile,xpolcoffset,"UnsignedInteger16",ByteOrdering->+1];
BinaryWrite[outfile,dataCadj[[j,i]],"UnsignedInteger16",ByteOrdering->+1];
}
,{i,1,chartchannels}
];
BinaryWrite[outfile,65278,"UnsignedInteger16",ByteOrdering->+1];
}
,{j,1,chartrecs}
];
Close[outfile];
Print["Xpol spectrograph data file saved as "<>outfilename]
 *)

(* 


Do data decimation FIRST to keep array size small as possible.

Instead of making raster arrays containing the color triplets,
Try ArrayPlot[] with ColorFunction-->ColorsInT[[data point value]]
eg. ColorFunction->Function[{x,y,z},Hue[z]]
ColorFunction->Function[{a},ColorsInT[[a]]]
The use of ArrayPlot[] will allow the use of tick marks, etc to do axis labeling 
or
ColorFunction->Function[{a},RGBColor[ColorsInT[[a,1]],ColorsIn[[a,2]],ColorsIn[[a,3]]]]


*)




output=Table["",{i,1,500}];
output[[1]]=Style["PDA (Period Domain Algorithm):"<>FromCharacterCode[13]<>"A Method of Searching for Pulsar Signals in Radio-SkyPipe Data",Larger,Darker[Darker[Blue]],Bold];
output[[2]]="Dave Typinski,2009";
output[[4]]="Job started "<>ToString[CurrentDate]<>", "<>ToString[CurrentTime];
output[[10]]=Style["PDA Analysis",Bold];
output[[12]]="                   Input File name: "<>filename;
output[[14]]="     Search period center (P0) (s): "<>ToString[PaddedForm[centerperiod,{15,13},NumberPadding->{"","0"}]];
output[[15]]="           Search period step (μs): "<>ToString[deltat*10^6];
output[[16]]="                 Number of periods: "<>ToString[numpers];
output[[17]]="                 Search range (μs): "<>ToString[periodbandwidth-deltat*10^6];
output[[18]]="           Search period start (s):"<>ToString[PaddedForm[periodstart,{8,7}]];
output[[19]]="             Search period end (s):"<>ToString[PaddedForm[periodend,{8,7}]];
output[[20]]="                    Number of bins: "<>ToString[bins];
output[[22]]="              PDA mean (K x 1,000): "<>ToString[PaddedForm[Round[Mean[averagebinsstatsplot[[All,2]]]*10^3]/10.^3,{5,3},NumberPadding->{"","0"}]];
charttitle=Column[{Style["PDA Analysis",Plain,FontSize->12,FontFamily->Helvetica],
Style["File: "<>filename,Plain,FontSize->11,FontFamily->Helvetica]},Center];
plot1=ListPlot[averagebinsstatsplot,
PlotRange->Full,
PlotStyle->{RGBColor[0,0,.6],Thin},
PlotMarkers->{Style["+",Bold,FontSize->22,FontFamily->Helvetica]},
Axes->None,
GridLines->None,
Frame->True,
RotateLabel-> True,FrameStyle->{{Directive[Thin,FontFamily->Helvetica,11],Thin},{Directive[Thin,FontFamily->Helvetica,10],Thin}},FrameLabel->{Style["Period Offset (μs) from "<>ToString[PaddedForm[centerperiod,{15,13},NumberPadding->{"","0"}]]<>" s",FontSize->10,FontFamily->Helvetica],Style["Maximum Output Bin Temperature (K x 1,000)",FontSize->11,FontFamily->Helvetica],charttitle,""}];
output[[99]]=Rasterize[plot1,ImageSize->650];
(* Export raster output to PDF in same directory as data file *)
g=Panel[Column[{Column[output[[10;;11]]],Column[output[[12;;31]],Alignment->":"],output[[99]]}],Column[output[[1;;5]]]];
m=Panel[output[[199]]];
h=Panel[Column[{Column[output[[201;;202]]],Column[output[[203;;212]],Alignment->":"],output[[213]]}]];
i=Panel[Column[{Column[output[[301;;302]]],Column[output[[303;;309]],Alignment->":"]}]];
(*j=Column[SFAout[[All]]];*)
k=Column[{g,m,h,i}];
outfname=filepath<>"PDA Bins="<>ToString[bins]<>" P0="<>ToString[PaddedForm[centerperiod,{7,6},NumberPadding->{"","0"}]]<>" step="<>ToString[deltat*10^6]<>"us for "<>StringDrop[filename,-3]<>"pdf";
 Export[outfname,k];
g
m
h
i