Commit 9ee1b0c5 authored by Timothy Shippert's avatar Timothy Shippert
Browse files

Update to langley to work with mfrsr7nch

* Process langley_mfrsr7nch
* Lots of switching on whether broadband field exists or not
* Some refactoring of channel ordering and number of narrowband channels
* Refactor plotting routine to work with  7nch and also plot actual_wavelength

Passes current langley release tests for bb/6nch configuration.  Runs for
7nch dataset (2018-2020), but output not verified yet.
parent 98b69343
......@@ -114,10 +114,10 @@ DATA* barnard_langley (DATA *D, float *gNomCal)
* we will extract the "Input_Datastreams" information
*/
for (i=0; i < D->nSamples[B1][o][0][IN_BROADBAND]; i++, ot++)
for (i=0; i < D->nSamples[B1][o][0][IN_CHN1]; i++, ot++)
{
if (D->sampleTimes[B1][o][0][IN_BROADBAND][i].zt_Sec >= startime &&
D->sampleTimes[B1][o][0][IN_BROADBAND][i].zt_Sec <= etime)
if (D->sampleTimes[B1][o][0][IN_CHN1][i].zt_Sec >= startime &&
D->sampleTimes[B1][o][0][IN_CHN1][i].zt_Sec <= etime)
{
if (sot == -1)
{
......@@ -156,14 +156,14 @@ DATA* barnard_langley (DATA *D, float *gNomCal)
ot = 0;
for (o=0; o < D->nObs[B1]; o++)
{
for (i=0; i < D->nSamples[B1][o][0][IN_BROADBAND]; i++, ot++)
for (i=0; i < D->nSamples[B1][o][0][IN_CHN1]; i++, ot++)
{
if (ot >= sot && ot <= eot)
{
/* do the time first - we need to calculate */
/* efactor from first sample time anyway */
/* sample time */
ztime[n] = D->sampleTimes[B1][o][0][IN_BROADBAND][i];
ztime[n] = D->sampleTimes[B1][o][0][IN_CHN1][i];
TC_ZtToUI(&ztime[n], &t);
jday = yymmdd_to_juldate(t.ds_yymmdd);
year = t.ds_yymmdd / 10000;
......@@ -251,7 +251,7 @@ DATA* barnard_langley (DATA *D, float *gNomCal)
} /* end else if (A[n] >= ... */
else if ((A[n] <= HIGH_AM && A[n] >= LOW_AM) &&
(o == (D->nObs[B1] - 1) &&
i == D->nSamples[B1][o][0][IN_BROADBAND] - 1))
i == D->nSamples[B1][o][0][IN_CHN1] - 1))
{
/*************************************************************/
/* this means that our last sample was not outside the valid */
......@@ -287,17 +287,31 @@ DATA* barnard_langley (DATA *D, float *gNomCal)
}
} /* if (n > ... */
/////////////////////////////////////////////////////////
/* Finally, the ln(radiance/efactor) */
/* broadband channel in seperate field */
lnI[0][n] = D->BWdata[B1][o][0][IN_BROADBAND][i][0] > 0 ?
log (D->BWdata[B1][o][0][IN_BROADBAND][i][0] /
efactor) : -9999;
// TRS 20191230
// If broadband exists, put it first, otherwise, just loop
// over NCHANN
int kstart=1; // first narrowband channel in lnI[k] field
if (D->nSamples[B1][o][0][IN_BROADBAND] == 0) {
// This means we have 7 narrowband channels
kstart=0;
} else {
// This means we have a broadband channel in seperate
// field, and tack it onto the front
lnI[0][n] = D->BWdata[B1][o][0][IN_BROADBAND][i][0] > 0 ?
log (D->BWdata[B1][o][0][IN_BROADBAND][i][0] /
efactor) : -9999;
}
/* narrowband channels */
for (k=1; k < NCHANNELS; k++)
for (k=kstart; k < NCHANNELS; k++)
{
lnI[k][n] = D->BWdata[B1][o][0][IN_CHN1 + (k-1)][i][0] > 0.0 ?
log (D->BWdata[B1][o][0][IN_CHN1 + (k-1)][i][0] /
lnI[k][n] = D->BWdata[B1][o][0][IN_CHN1 + (k-kstart)][i][0] > 0.0 ?
log (D->BWdata[B1][o][0][IN_CHN1 + (k-kstart)][i][0] /
efactor) : -9999;
} /* end for (k... */
......
......@@ -234,7 +234,11 @@ long etime;
char gSite[20];
char gFacility[20];
void getWavelengthAttr ();
// This is used to toggle broadband vs. filter 7 configuration
// We make it global so we can use it in create_quicklooks
int broadband_exists=0;
void getWavelengthAttr (DATA *, int);
DATA* ComputeNewData (DATA *D)
/*****************************************************************************
......@@ -261,7 +265,7 @@ DATA* ComputeNewData (DATA *D)
int getStartEndTime3 (DATA *, long *, long *);
DATA *barnard_langley (DATA *, float *);
int michalsky_langley (DATA *, DATA *, float *);
int writeWavelengthAttr (DATA *);
int writeWavelengthAttr (DATA *, int);
char *head_id;
char *logger_id;
......@@ -271,7 +275,6 @@ DATA* ComputeNewData (DATA *D)
void checkQCFlag();
float gNomCal[7];
/* First thing, store global variables gSite and gFacility */
sprintf (gSite, "%s", dsproc_get_site());
......@@ -284,19 +287,33 @@ DATA* ComputeNewData (DATA *D)
/* We did that in parse_commandline_options ... so, we will always want to use "p = 0" */
B1 = 0;
// Figure out if we have a broadband channel or a filter7 channel
if (D->nSamples[B1][0][0][IN_BROADBAND] > 0) {
broadband_exists=1;
} else {
broadband_exists=0;
}
/* Convert input data back to counts */
for (o=0; o < D->nObs[0]; o++)
{
/* p o [0] f t n */
memcpy(&gNomCal[0], (float *) &D->BWdata[B1][o][0][14][0][0], sizeof(float));
memcpy(&gNomCal[1], (float *) &D->BWdata[B1][o][0][15][0][0], sizeof(float));
memcpy(&gNomCal[2], (float *) &D->BWdata[B1][o][0][16][0][0], sizeof(float));
memcpy(&gNomCal[3], (float *) &D->BWdata[B1][o][0][17][0][0], sizeof(float));
memcpy(&gNomCal[4], (float *) &D->BWdata[B1][o][0][18][0][0], sizeof(float));
memcpy(&gNomCal[5], (float *) &D->BWdata[B1][o][0][19][0][0], sizeof(float));
memcpy(&gNomCal[6], (float *) &D->BWdata[B1][o][0][20][0][0], sizeof(float));
// This is not the right way to do this, but we'll assume, for now
// that the calibration factors fields come right after the QC
// fields. If we have a broadband channel, then start there,
// otherwise skip one to start at filter 1.
int CALIBx = QC_CHN7+1;
if (! broadband_exists) CALIBx += 1;
memcpy(&gNomCal[0], (float *) &D->BWdata[B1][o][0][CALIBx][0][0], sizeof(float));
memcpy(&gNomCal[1], (float *) &D->BWdata[B1][o][0][CALIBx+1][0][0], sizeof(float));
memcpy(&gNomCal[2], (float *) &D->BWdata[B1][o][0][CALIBx+2][0][0], sizeof(float));
memcpy(&gNomCal[3], (float *) &D->BWdata[B1][o][0][CALIBx+3][0][0], sizeof(float));
memcpy(&gNomCal[4], (float *) &D->BWdata[B1][o][0][CALIBx+4][0][0], sizeof(float));
memcpy(&gNomCal[5], (float *) &D->BWdata[B1][o][0][CALIBx+5][0][0], sizeof(float));
memcpy(&gNomCal[6], (float *) &D->BWdata[B1][o][0][CALIBx+6][0][0], sizeof(float));
}
/* assume all obs are at the same location */
......@@ -324,7 +341,7 @@ DATA* ComputeNewData (DATA *D)
NULL,
NULL);
getWavelengthAttr (D);
getWavelengthAttr (D, broadband_exists);
/* bbb */
......@@ -362,7 +379,7 @@ DATA* ComputeNewData (DATA *D)
NULL,
NULL);
getWavelengthAttr (D);
getWavelengthAttr (D, broadband_exists);
if ((newD = barnard_langley (D, (float *) &gNomCal)) == NULL) {
bw_exit(newD, "Problem in Barnard; continuing to next processing interval");
}
......@@ -371,7 +388,7 @@ DATA* ComputeNewData (DATA *D)
bw_exit(newD, "Problem in Michalsky; continuing to next processing interval");
};
writeWavelengthAttr (newD);
writeWavelengthAttr (newD, broadband_exists);
for (o=0; o < newD->nObs[LANGLEY]; o++)
{
......@@ -552,6 +569,7 @@ void create_quicklooks (void)
char command[1024], argument[1024];
char setdate[1024],
setplatform[1024], setsite[1024], setfacility[1024];
char buf[1024];
char **argv, *vapname;
int argc;
......@@ -584,6 +602,10 @@ void create_quicklooks (void)
* Set some environment variables to guide IDL...
*/
// First of all, do we have a broadband channel or filter7?
sprintf(buf,"BROADBAND=%d", broadband_exists);
putenv(buf);
/*
* It is important that each environment variable has its own
* memory allocated to it (ie., its own variable). This has
......@@ -596,13 +618,23 @@ void create_quicklooks (void)
/* For the ADI version of this VAP, we need to pull out either
* "nimfr" or "mfrsr" from the platform
*
* TRS 20200106 - need to check for 7nch first, mfrsr will match
* the mfrsr7nch as well
*
*/
if (strstr(platform, "mfrsr") != (char *) NULL)
if (strstr(platform, "mfrsr7nch") != (char *) NULL)
{
sprintf (setplatform, "PLATFORM=mfrsr7nch");
}
else if (strstr(platform, "mfrsr") != (char *) NULL)
{
sprintf (setplatform, "PLATFORM=mfrsr");
}
else if (strstr(platform, "nimfr7nch") != (char *) NULL)
{
sprintf (setplatform, "PLATFORM=nimfr7nch");
}
else
{
sprintf (setplatform, "PLATFORM=nimfr");
......@@ -888,16 +920,26 @@ void checkQCFlag (DATA *D)
} /* end checkQCFlag (... */
void getWavelengthAttr (DATA *D)
void getWavelengthAttr (DATA *D, int broadband_exists)
{
char *foo;
int r;
int r, rb;
// rn is number of broadband channels, to subtract from total number of
// channels when indexing.
if (broadband_exists) {
rb=1;
} else {
rb=0;
}
/*
* it is assumed that the chns will not change over obs.
* if they do then this does not work and will need to be rewritten
*/
for (r = 0; r < (NCHANNELS - 1); r++)
// We have NCHANNELS-rb narrowband channels
for (r = 0; r < NCHANNELS-rb; r++)
{
foo = (char *)bw_GetAttrValue (D,
B1,
......@@ -919,38 +961,40 @@ void getWavelengthAttr (DATA *D)
if (foo != NULL)
{
channels[r + 1] = atof (foo);
// Note that we shift one up to hold broadband, if it exists
channels[r + rb] = atof (foo);
// printf ("DEBUG: r = %d ... actual_wavelength = %f\n", r, atof(foo));
}
} /* end for (r... */
channels[0] = 0; /* This is the broadband channel */
if (broadband_exists) {
channels[0] = 0; /* Broadband channel in front */
}
/* return (TRUE); */
} /* end getWavelengthAttr (... */
int writeWavelengthAttr (DATA *newD)
int writeWavelengthAttr (DATA *newD, int broadband_exists)
{
char chns[16];
int o, r;
int o, r, r0;
if (broadband_exists) {
// Broadband mfrsr config
// Don't write out actual_wavelength to broadband channel (r=0)
r0=1;
} else {
// 7nch config, so write out all seven actual_wavelengths
r0=0;
}
for (o = 0; o < newD->nObs[LANGLEY]; o++)
{
for (r = 0; r < NCHANNELS; r++)
for (r = r0; r < NCHANNELS; r++)
{
/* the bb chn is the last chn in the array */
if (r == 0)
{
continue; /* Don't set the actual_wavelength value for the
* broadband channel
*/
}
else
{
sprintf (chns, "%.3f nm", channels[r - 1]);
}
sprintf (chns, "%.3f nm", channels[r - r0]);
/* barnard */
bw_SetAttrValue (newD, LANGLEY, o, OD + r, DCT_String,
......@@ -988,19 +1032,9 @@ int writeWavelengthAttr (DATA *newD)
for (o = 0; o < newD->nObs[PLOT]; o++)
{
for (r = 0; r < NCHANNELS; r++)
for (r = r0; r < NCHANNELS; r++)
{
/* the bb chn is the last chn in the array */
if (r == 0)
{
continue; /* Don't set the actual_wavelength value for
* the broadband channel
*/
}
else
{
sprintf (chns, "%f nm", channels[r - 1]);
}
sprintf (chns, "%f nm", channels[r - r0]);
/* barnard */
bw_SetAttrValue (newD, PLOT, o, LN_I + r, DCT_String,
......
......@@ -69,6 +69,11 @@ int B1;
#define NPRESPLATS 4
#define PRES_CONF_FILE "langley_pres.conf"
// TRS 20191230 - update this to include both broadband and chn7, because
// we need langley_retriever.h to include both, so it will work on both
// legacy mfrsr and nch7 mfrsr. This gives us the indices - the BWdata
// structure for missing input fields should just be a null pointer, but
// this way the indices stay the same no matter what.
enum input_fields
{
IN_BROADBAND = 0,
......@@ -78,13 +83,15 @@ int B1;
IN_CHN4,
IN_CHN5,
IN_CHN6,
IN_CHN7,
QC_BROADBAND,
QC_CHN1,
QC_CHN2,
QC_CHN3,
QC_CHN4,
QC_CHN5,
QC_CHN6
QC_CHN6,
QC_CHN7
}; /* enum input_fields */
enum input_pres_flds
......
......@@ -71,6 +71,10 @@
;
; quiet - If set, don't print status updates or names of output
; files; will still print final stats if /stats is set.
;
; broadband - if set to 1, we are using original broadband channel
; configuration; if 0, we are using 7nch configuration
; (filter7 at 1624 nm)
;
; plot_to - sets the type of output plots. Can be set to ps (postscript),
; png (png graphics file), or x (to the X display). Default is
......@@ -99,7 +103,7 @@ pro langley, ch, to = to, date = td, $
algorithm = algorithm, uncalibrated = uncalibrated, $
onefile = onefile, noplot = noplot, stats = stats, $
quiet = quiet, plot_to = plot_to, no_catch = no_catch, $
zlog = zlog
zlog = zlog, broadband = broadband
;
; force this rascal to use the Z buffer (so quicklooks will
......@@ -190,7 +194,12 @@ pro langley, ch, to = to, date = td, $
return
endif
endif
endif
;; Broadband on by default
if (not keyword_set(broadband)) then begin
broadband=1
endif
; Revision number of this file. It shows up in tiny print on the Langley
; plots themselves.
......@@ -308,9 +317,15 @@ pro langley, ch, to = to, date = td, $
; These shouldn't be hardcoded; we should be able to pull it out from
; the data itself. Oh well.
; Should move this down to the channel loop, and extract actual_wavelength
NCHNS = 7
if (broadband ne 1) then begin
NCHNS = 8 ; not really, but it's just easier if we spoof index 0 as broadband
endif
chname = ['broadband', '415 nm', '499 nm', '608 nm', '664 nm', $
'860 nm', '938 nm']
'860 nm', '938 nm', '1624 nm']
if (algor eq MICHALSKY) then begin
rejwhy = ['no points to regress', 'all points filtered', $
......@@ -342,6 +357,9 @@ pro langley, ch, to = to, date = td, $
if (ch le 0) then begin
ch = 0
endch = 6
if (broadband eq 0) then begin
endch = 7
endif
endif else begin
endch = ch
endelse
......@@ -392,52 +410,61 @@ pro langley, ch, to = to, date = td, $
ncdf_varget, fin, algorithm + '_airmass', airmass
if (algorithm eq 'barnard') then begin
ncdf_varget, fin, 'barnard_lnI_broadband', foonew
endif
if (algorithm eq 'michalsky') then begin
ncdf_varget, fin, 'michalsky_lnI_broadband', foonew
endif
lograd[*,0] = foonew
foonew = 0
if (algorithm eq 'barnard') then begin
ncdf_varget, fin, 'barnard_rejected_broadband', foonew
rejected[*,0] = foonew
endif
foonew = 0
if (broadband eq 1) then begin
chname[0] = 'broadband'
if (algorithm eq 'barnard') then begin
ncdf_varget, fin, 'barnard_lnI_broadband', foonew
endif
if (algorithm eq 'michalsky') then begin
ncdf_varget, fin, 'michalsky_lnI_broadband', foonew
endif
lograd[*,0] = foonew
foonew = 0
if (algorithm eq 'barnard') then begin
ncdf_varget, fin, 'barnard_rejected_broadband', foonew
rejected[*,0] = foonew
endif
foonew = 0
if (algorithm eq 'michalsky') then begin
ncdf_varget, fin, 'michalsky_rejected_broadband', foonew
rejected[*,0] = foonew
endif
foonew = 0
if (algorithm eq 'michalsky') then begin
ncdf_varget, fin, 'michalsky_rejected_broadband', foonew
rejected[*,0] = foonew
endif
foonew = 0
for j=1, NCHNS - 1 do begin
if (algorithm eq 'barnard') then begin
ncdf_varget, fin, 'barnard_lnI_filter' + strtrim (j, 2), foonew
endif
if (algorithm eq 'michalsky') then begin
ncdf_varget, fin, 'michalsky_lnI_filter' + strtrim (j, 2), foonew
endif
field = algorithm + '_lnI_filter' + strtrim(j,2)
ncdf_varget, fin, field, foonew
lograd[*, j] = foonew
foonew = 0
;; Want to set the chname to use the actual wavelength
;; GRR. Verifying a field attribute exists is so hard
;; that I'm going to do it anyway, just out of spite
ai=ncdf_varinq(fin, field)
for a=0,ai.natts-1 do begin
aname=ncdf_attname(fin, field, a)
if (aname eq 'actual_wavelength') then begin
ncdf_attget,fin,field,'actual_wavelength',awl
;; Seriously - byte to string to float to string, to
;; format it correctly. IDL is dumn
chname[j] = strtrim(string(float(string(awl)),f='%6.1f'),2)
endif
endfor
if (algorithm eq 'barnard') then begin
ncdf_varget, fin, 'barnard_rejected_filter' + strtrim (j, 2), foonew
endif
if (algorithm eq 'michalsky') then begin
ncdf_varget, fin, 'michalsky_rejected_filter' + strtrim (j, 2), foonew
endif
foonew = 0
field = algorithm + '_rejected_filter' + strtrim(j,2)
ncdf_varget, fin, field, foonew
rejected[*, j] = foonew
foonew = 0
......@@ -515,32 +542,37 @@ pro langley, ch, to = to, date = td, $
endif else begin
scvar = 'barnard_solar_constant_sdist'
sdvar = 'barnard_error_fit'
ncdf_varget, flang, 'barnard_error_slope_broadband', foo
sigs[*,0] = foo
foo = 0
if (broadband eq 1) then begin
ncdf_varget, flang, 'barnard_error_slope_broadband', foo
sigs[*,0] = foo
foo = 0
endif
endelse
ncdf_varget, flang, algorithm + '_optical_depth_broadband', foo
tau[*, 0] = foo
foo = 0
ncdf_varget, flang, scvar + '_broadband', foo
I0[*, 0] = foo
foo = 0
ncdf_varget, flang, sdvar + '_broadband', foo
sd[*, 0] = foo
foo = 0
ncdf_varget, flang, algorithm + '_number_points_broadband', foo
npoints[*,0] = foo
foo = 0
ncdf_varget, flang, algorithm + '_good_fraction_broadband', foo
frac[*, 0] = foo
foo = 0
ncdf_varget, flang, algorithm + '_badflag_broadband', foo
bad[*, 0] = foo
foo = 0
if (broadband eq 1) then begin
ncdf_varget, flang, algorithm + '_optical_depth_broadband', foo
tau[*, 0] = foo
foo = 0
ncdf_varget, flang, scvar + '_broadband', foo
I0[*, 0] = foo
foo = 0
ncdf_varget, flang, sdvar + '_broadband', foo
sd[*, 0] = foo
foo = 0
ncdf_varget, flang, algorithm + '_number_points_broadband', foo
npoints[*,0] = foo
foo = 0
ncdf_varget, flang, algorithm + '_good_fraction_broadband', foo
frac[*, 0] = foo
foo = 0
ncdf_varget, flang, algorithm + '_badflag_broadband', foo
bad[*, 0] = foo
foo = 0
endif
for j=1, NCHNS - 1 do begin
ncdf_varget, flang, algorithm + '_optical_depth_filter' + $
strtrim (j, 2), foo
strtrim (j, 2), foo
tau[*, j] = foo
foo = 0
ncdf_varget, flang, scvar + '_filter' + strtrim (j, 2), foo
......
......@@ -21,15 +21,20 @@ pro langley_batch
site = getenv ("SITE")
facility = getenv ("FACILITY")
algorithm = getenv ("ALGORITHM")
broadband = getenv ("BROADBAND")
resolve_routine, 'langley'
; resolve_all
langley, 1, to=6, date = date, platform = platform, site=site, $
facility=facility, algorithm = 'barnard', $
; Broadband, plot channels 1-6, otherwise plot to 7
bb=fix(broadband)
langley, 1, to=7-bb, date = date, platform = platform, site=site, $
facility=facility, algorithm = 'barnard', broadband = broadband, $
plot_to = 'png'
langley, 1, to=6, date = date, platform = platform, site=site, $
facility=facility, algorithm = 'michalsky', $
langley, 1, to=7-bb, date = date, platform = platform, site=site, $
facility=facility, algorithm = 'michalsky', broadband = broadband, $
plot_to = 'png'
; langley, 1, to=6, date = date, platform = platform, site=site, $
; facility=facility, algorithm = 'barnard', /onefile, $
......
......@@ -23,6 +23,7 @@ const char *Field_List[][BW_MAX_FIELDS]=
"direct_normal_narrowband_filter4",
"direct_normal_narrowband_filter5",
"direct_normal_narrowband_filter6",
"direct_normal_narrowband_filter7",
"qc_direct_normal_broadband",
"qc_direct_normal_narrowband_filter1",
"qc_direct_normal_narrowband_filter2",
......@@ -30,6 +31,7 @@ const char *Field_List[][BW_MAX_FIELDS]=
"qc_direct_normal_narrowband_filter4",
"qc_direct_normal_narrowband_filter5",
"qc_direct_normal_narrowband_filter6",
"qc_direct_normal_narrowband_filter7",
"nominal_calibration_factor_broadband",
"nominal_calibration_factor_filter1",
"nominal_calibration_factor_filter2",
......@@ -37,6 +39,7 @@ const char *Field_List[][BW_MAX_FIELDS]=
"nominal_calibration_factor_filter4",
"nominal_calibration_factor_filter5",
"nominal_calibration_factor_filter6",
"nominal_calibration_factor_filter7",
NULL
},
};
......
......@@ -99,10 +99,10 @@ int setup_for_michalsky (DATA *D)
sot = -1;
for (o=0; o<D->nObs[B1]; o++)