Skip to content

GitLab

Commit 02287497 authored by Chitra Sivaraman's avatar Chitra Sivaraman
Browse files

Initial files to code.arm.gov.

parents
Hide whitespace changes
Inline Side-by-side
Showing with 7168 additions and 0 deletions
src/barnard_langley.c 0 → 100644
View file @ 02287497
/*******************************************************************************
* COPYRIGHT (C) 2001 Battelle Memorial Institute.
* All Rights Reserved. (Now, with this out of the way...)
*
* RCS INFORMATION:
* $RCSfile: barnard_langley.c,v $
* $Revision: 1.19 $
* $Author: koontz $
* $Locker: $
* $Date: 2011/12/29 18:56:44 $
* $State: Exp $
* $Name: $
* $Id: barnard_langley.c,v 1.19 2011/12/29 18:56:44 koontz Exp $
*
* FUNCTIONS IN THIS FILE:
* barnard_langely (DATA *, DATA *)
*
* DESIGN:
* <Discuss overall design/purpose of class.>
*******************************************************************************/
#define BARNARD_LANGLEY_C
/****** General Includes ******/
#include <stdio.h>
#include <string.h>
#include <math.h>
/****** Zebra Includes ******/
/* #include "defs.h"
* #include "message.h"
* #include "DataStore.h"
*/
/****** BW Includes ******/
/* #define BW_CODE
* #include "bw_main.h"
* #include "bw_retriever.h"
*/
#define BW_CODE
#include "bw_adi.h"
/****** Application Includes ******/
#include "langley.h"
/*******************************************************************************
* Function:
* barnard_langley -
*
* Inputs:
*******************************************************************************/
DATA* barnard_langley (DATA *D, float *gNomCal)
{
DATA *newD;
ZebTime ztime[4325];
UItime t;
int ot, sot, eot;
int f, i,j,k, l, m, n,o, s, nsamples, npoints, norig, nruns, nmax;
int nrej, nzero, old_rej;
int nobs[NOUTPLATS] = {1, 2};
int badflag;
int hh, mm, ss, jday, old_jday;
int start[NSLABS], end[NSLABS];
float frac;
float lnI[NCHANNELS][4325], A[4325], AZ[4325],
lnItemp[4325], Atemp[4325],
weights[4325], wtemp[4325];
double airmass();
/* time stuff */
int year;
double dayfrac, arg, az, distance;
float sdist[4325];
float sdist_avg;
int sdist_count;
double efactor, el, ap_ra, ap_dec, refrac;
extern int solarposition();
/* fit variables */
float sigp, intercept, slope, sigi, sigs, y;
float resid;
#ifdef OLDFIT
extern fit_ ();
#endif
void lfit();
int yymmdd_to_juldate();
extern long startime;
extern long etime;
float dvalue;
/* function prototypes */
DATA* maker_AllocateMemoryForDATA();
nsamples = 0;
efactor = -9999; /* initialize */
/* determine the number of samples */
ot = 0;
eot = 0;
sot = -1;
for (o=0; o < D->nObs[B1]; o++)
{
/* Get the "ingest_software" info for this observation, from which
* we will extract the "Input_Datastreams" information
*/
for (i=0; i < D->nSamples[B1][o][0][IN_BROADBAND]; 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 (sot == -1)
{
sot = ot;
/* Using this observation.
*
* Get version from a global attribute in the observation
* of MFR data we are using for the Langley.
*/
}
eot = ot;
} /* end if (D->sampleTimes... */
} /* end for (t... */
} /* end for (o... */
if (sot == -1)
{
bw_return(NULL, "Barnard: no samples for given time period");
}
nsamples = (eot - sot) + 1;
/******************************************************/
/* Now pull the BW fields into our lnI arrays - we */
/* have to multiply by efactor and take the log */
/* anyway, so we might as well loop over all obs and */
/* make one big array of all our values. Airmasses, */
/* too. */
/******************************************************/
n = 0;
m = 0;
old_jday = -9999;
ot = 0;
for (o=0; o < D->nObs[B1]; o++)
{
for (i=0; i < D->nSamples[B1][o][0][IN_BROADBAND]; 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];
TC_ZtToUI(&ztime[n], &t);
jday = yymmdd_to_juldate(t.ds_yymmdd);
year = t.ds_yymmdd / 10000;
if (year < 1900) year += 1900; /* y2k or not */
hh = t.ds_hhmmss / 10000;
mm = (t.ds_hhmmss /100 ) % 100;
ss = t.ds_hhmmss % 100;
dayfrac = (double) (hh*3600 + mm*60 + ss) / (24.0*3600.0);
dayfrac += (double) jday;
/* get efactor, if it's a new day */
if (efactor <= -8888 || jday != old_jday)
{
/* get efactor from first sample */
if (year % 4 == 0)
arg = (double) (jday - 1)/366.0*2.0*M_PI;
else
arg = (double) (jday - 1)/365.0*2.0*M_PI;
efactor = 1.00011 + .034221*cos(arg) + .00128*sin(arg) +
.000719*cos(2*arg) + 7.7e-5*sin(2*arg);
} /* end if (efactor... */
/* now get airmasses */
solarposition(year, 0, dayfrac, (double) 0.0, lat, lon,
&ap_ra, &ap_dec, &el, &refrac, &az, &distance);
sdist[n] = distance;
A[n] = airmass(el+refrac);
AZ[n] = az;
/* don't do it for first point... */
if (n > 1)
{
/* weight by inverse density of points, that is, by the */
/* gradiant of airmass vs. time */
/* we need the n+1 airmass and time to do the nth weight */
/* so we fill in the j=n-1 weight on the nth iteration */
j = n - 1;
#define TWEIGHT
#ifdef TWEIGHT
/* a scale factor of 2000.0 sets the weights in the middle*/
/* of a typical airmass slab to just about 1. This ought */
/* to help deal with our precision issues. */
// Check to see if our neighboring airmasses are bad -
// this should only happen when there's a gap in the
// data, as we only fit for A in [2,6]. Anyway, I'll use
// < 0 and > 10*HIGH_AM as the cutoff
if (A[j+1] < 0 || A[j+1] > 10*HIGH_AM ||
A[j-1] < 0 || A[j-1] > 10*HIGH_AM) {
weights[j]=0.0;
} else {
weights[j] = 2000.0*fabs((A[j+1] - A[j-1])/
(float) (ztime[j+1].zt_Sec - ztime[j-1].zt_Sec));
}
#else
weights[j] = 1.0;
#endif
/* mark start and end of our airmass slabs within range */
if (A[n] >= LOW_AM && A[n] <= HIGH_AM &&
(A[n-1] < LOW_AM || A[n-1] > HIGH_AM))
{
/* index start of slab */
start[m] = n;
} /* end if (A[n] >= ... */
else if (A[n] >= LOW_AM && A[n] <= HIGH_AM &&
(AZ[n] >= 180.0 && AZ[n - 1] < 180.0))
{
/* index end and start at solar noon local time */
end[m++] = n - 1;
start[m] = n;
} /* end else if *AZ[n] >= ... */
else if ((A[n-1] >= LOW_AM && A[n-1] <= HIGH_AM &&
(A[n] < LOW_AM || A[n] > HIGH_AM)))
{
/* index end of slab */
end[m++] = n - 1;
} /* 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))
{
/*************************************************************/
/* this means that our last sample was not outside the valid */
/* airmass range - in other words, we probably didn't get all*/
/* the data for this day. We'll go ahead and peg the end of */
/* the slab and do a langley plot anyway - but we send a */
/* warning message, too. */
/*************************************************************/
msg_ELog(EF_PROBLEM,
"==> MISSING some daylight data for final Langley plot.");
msg_ELog(EF_PROBLEM,
"==> Possibly next day's data is missing.");
msg_ELog(EF_PROBLEM,
"==> May have to rerun when data is available to fix.");
end[m++] = n;
} /* end else if ((A[n] <= ... */
else if ((A[n] <= HIGH_AM && A[n] >= LOW_AM) && ot==eot)
{
// This means we are still in the Airmass range but
// we've hit the end of our analysis period. If the
// above test didn't catch it, it probably means that
// eot is local midnight, but the sun is still up,
// which probably means arctic or antarctic. At any
// rate, we have to close our langley slab, so:
msg_ELog(EF_INFO,
"==> Last airmass in slab %f still in range [%f, %f]",
A[n], LOW_AM, HIGH_AM);
msg_ELog(EF_INFO,
"==> Probably arctic or antarctic and sun is up at midnight");
msg_ELog(EF_INFO,
"==> Continuing with Langley analysis for this slab");
end[m++] = n;
}
} /* 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;
/* narrowband channels */
for (k=1; 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] /
efactor) : -9999;
} /* end for (k... */
/* increment lnI counter */
old_jday = jday;
n++;
} /* end if (ot... */
} /* end for (i... */
} /* end for (o... */
/* setup our newD to accept data */
nobs[LANGLEY] = 1; /* writes out one file per run */
#ifdef PLOTDATA
nobs[PLOT] = m; /* one plotfile for each slab */
if (m <= 0)
bw_return(NULL, "==> No usable data for this time window (probably due to lack of daylight) <==");
#endif
nmax = 0;
for (i=0; i<m; i++)
{
npoints = end[i] - start[i] + 1;
if (npoints > nmax) nmax = npoints;
} /* end for (i... */
newD = (DATA *) maker_AllocateMemoryForDATA (nobs, nmax);
/***********************************************/
/* Okay, now we have big long lnI and A arrays */
/* Now we loop over each airmass slab within */
/* [2,6] and do another set of fits */
/***********************************************/
printf ("Doing main Barnard storage loop ... m= %d\n", m);
sdist_avg = 0.0;
sdist_count = 0;
for (i = 0; i < m; i++)
{
/* number of points in this slab */
npoints = end[i] - start[i] + 1;
norig = npoints;
nrej = 0;
for (k=start[i]; k < end[i]; k++)
{
sdist_avg = sdist_avg + sdist[k];
sdist_count++;
}
/* need two temp arrays, so we can shuffle back and forth */
for (k=0,l=0;k<NCHANNELS;k++)
{
/* initialize the temp arrays */
memcpy(lnItemp,&lnI[k][start[i]], npoints*sizeof(float));
memcpy(Atemp,&A[start[i]], npoints*sizeof(float));
memcpy(wtemp,&weights[start[i]], npoints*sizeof(float));
nrej = 0;
nzero=0;
/* scan down for missing radiances */
for (s=0;s<npoints;s++)
{
if (lnItemp[s] < -8888)
{
wtemp[s] = 0.0;
nzero++;
} /* end if (lnItemp[s]... */
} /* end for (s... */
#ifdef PLOTDATA
/* Write out our platform that contains the plot information */
for (s=0;s<npoints;s++)
{
/* things to do only once */
if (k == 0)
{
/* things to do for all fields */
for (f = 0; f < newD->nFields[PLOT][i]; f++)
{
/* write out sample times. */
newD->sampleTimes[PLOT][i][0][f][s] = ztime[start[i] + s];
/* write out nsamples */
newD->nSamples[PLOT][i][0][f] = npoints;
} /* end for (f... */
newD->BWdata[PLOT][i][0][AM][s][0] = Atemp[s];
/* write out sample location */
newD->subLoc[PLOT][i][0] = D->subLoc[B1][0][0];
/* write out new plot file for each slab */
newD->zebNewFileFlag[PLOT][i] = 1;
} /* end if (k... */
newD->BWdata[PLOT][i][0][LN_I+k][s][0] = lnItemp[s];
} /* end for (s... */
#endif /* end ifdef PLOTDATA */
/* now iterate over our fits */
nruns=0;
sigp = 1e9;
do
{
/* Here is our fit */
sigi=sigs=0;
lfit(Atemp, lnItemp, wtemp, npoints, &intercept, &slope,
&sigi, &sigs, &sigp, NULL);
#ifdef OLDFIT
fit_(Atemp, lnItemp, &npoints, sigtemp, &one, &intercept,
&slope, &sigi, &sigs, &chi2, &q);
sigp = sqrt(chi2/(float) (npoints-nrej-nzero-2));
/* we need to multiply sigs by sigp, because we've implicitly */
/* set our sigy=1 (in sigtemp) for all pts. in the fit. The */
/* estimated value for sigy is sigp... */
/* Really, I should just rewrite the fitting routine */
sigs *= sigp;
sigi *= sigp;
#endif /* end ifdef OLDFIT */
#define SIGLIM (k==1 ? .015 : .01)
old_rej = nrej;
/* Now run down the points and see if they are outside our fit */
for (j=0; j < npoints; j++)
{
y = intercept+slope*Atemp[j];
resid = lnItemp[j]-y;
/* reject data below the line and outside two sigmas */
if (wtemp[j] > 0 && resid < -2*sigp)
{
wtemp[j] = 0.0; /* this should zero out this point */
nrej++;
#ifdef PLOTDATA
/* tag the point as rejected in the plot platform */
newD->BWdata[PLOT][i][0][REJECTED+k][j][0] = 1;
#endif /* end ifdef PLOTDATA */
} /* end if (wtemp[j]... */
} /* end for (j... */
/* run the fit again, with the non rejected points */
nruns++;
} while (nruns < 20 && nrej > old_rej && sigp > SIGLIM);
/* take out the rejected points */
frac = (float)(npoints-nrej-nzero) / (float)npoints;
/* Now, set our badflags for various terms of rejection */
/* _track_ rejected fits */
if (sigp > SIGLIM)
badflag = 1;
else if (exp(intercept) < .05 ||
exp(intercept) > 5000.0)
badflag = 2;
else if (frac < .75 )
badflag = 3;
else
badflag = 0;
/* store our final field values */
if (isnan (slope) || isnan (intercept) || (npoints-nrej-nzero) < 3)
{
newD->BWdata[LANGLEY][0][0][OD + k][i][0] = MISSING;
newD->BWdata[LANGLEY][0][0][SC + k][i][0] = MISSING;
newD->BWdata[LANGLEY][0][0][ERR_FIT + k][i][0] = MISSING;
newD->BWdata[LANGLEY][0][0][ERR_SLOPE + k][i][0] = MISSING;
} /* end if (isnan (slope... */
else if (isinf ( slope) || isnan ( intercept) || (npoints-nrej-nzero) < 3)
{
newD->BWdata[LANGLEY][0][0][OD + k][i][0] = MISSING;
newD->BWdata[LANGLEY][0][0][SC + k][i][0] = MISSING;
newD->BWdata[LANGLEY][0][0][ERR_FIT + k][i][0] = MISSING;
newD->BWdata[LANGLEY][0][0][ERR_SLOPE + k][i][0] = MISSING;
} /* end if (isnan (slope... */
else if (badflag != 0)
{
/* The badflag value is not zero, so set fields missing */
newD->BWdata[LANGLEY][0][0][OD + k][i][0] = MISSING;
newD->BWdata[LANGLEY][0][0][SC + k][i][0] = MISSING;
newD->BWdata[LANGLEY][0][0][ERR_FIT + k][i][0] = MISSING;
newD->BWdata[LANGLEY][0][0][ERR_SLOPE + k][i][0] = MISSING;
}
else
{
newD->BWdata[LANGLEY][0][0][OD + k][i][0] = -slope;
dvalue = exp(intercept);
if (isnan(dvalue) || isinf(dvalue))
{
newD->BWdata[LANGLEY][0][0][SC + k][i][0] = MISSING;
}
else
{
/* For Filter 1 to 5 multiply the Solar Constant by gNomCal to convert the units
* back to "counts"
*/
dvalue = dvalue * gNomCal[k];
newD->BWdata[LANGLEY][0][0][SC + k][i][0] = dvalue;
}
/* final error in fit */
newD->BWdata[LANGLEY][0][0][ERR_FIT + k][i][0] = sigp;
/* error in slope */
newD->BWdata[LANGLEY][0][0][ERR_SLOPE + k][i][0] = sigs;
} /* end else */
newD->BWdata[LANGLEY][0][0][NPTS + k][i][0] = (float) (npoints-nrej-nzero);
newD->BWdata[LANGLEY][0][0][GOOD_FRACT + k][i][0] = frac;
newD->BWdata[LANGLEY][0][0][BAD_FLAG + k][i][0] = (float)badflag;
} /* end for (k... */
/* if key chn is good all chns are good */
if (newD->BWdata[LANGLEY][0][0][BAD_FLAG + KEY_CHN][i][0] == 0)
{
for (k = 0; k < NCHANNELS; k++)
newD->BWdata[LANGLEY][0][0][BAD_FLAG + k][i][0] = 0;
} /* end if (newD... */
else
{
for (k = 0; k < NCHANNELS; k++)
if (newD->BWdata[LANGLEY][0][0][BAD_FLAG + k][i][0] == 0)
newD->BWdata[LANGLEY][0][0][BAD_FLAG + k][i][0] = 4;
} /* end else */
/* report time in the middle of this slab. I do this way down */
/* here so I can trap on badflags in a later version if I want */
for (f=0; f<newD->nFields[LANGLEY][0]; f++)
{
newD->sampleTimes[LANGLEY][0][0][f][i]=ztime[start[i] + norig/2];
}
} /* end for (i... */
/***********************************************/
/* Okay, we now have values for all our slabs */
/* and for all channels. Just finish up with */
/* the boring stuff and we are ready to rock. */
/***********************************************/
/* Solar distance, Barnard ... may need to do some averaging ... */
dvalue = sdist_avg / sdist_count;
if (isinf(dvalue) || isnan(dvalue))
{
newD->BWdata[LANGLEY][0][0][BSDIST][0][0] = -9999.;
}
else
{
newD->BWdata[LANGLEY][0][0][BSDIST][0][0] = sdist_avg/sdist_count;
}
newD->nObs[LANGLEY] = 1;
for (f=0;f<newD->nFields[LANGLEY][0];f++)
{
newD->nSamples[LANGLEY][0][0][f] = m;
} /* end for (f... */
newD->subLoc[LANGLEY][0][0] = D->subLoc[B1][0][0];
/* platformNames = retriever_GetPlatformNames(); */
/* add a plot comment, if you want */
#if defined(TWEIGHT) && defined(PLOTDATA)
/* for (o=0; o<newD->nObs[PLOT]; o++) */
main_SetAttrValue (newD, PLOT, 0, -1, "plot_comment",
"Weighting by airmass gradient for fit");
#endif /* end if defined (... */
return (newD);
} /* barnard_langley */
int yymmdd_to_juldate (int yymmdd)
{
/* Local variables */
int i;
int leapyear;
int yy, mm, dd;
int remainder;
int month[12];
int jday = 0;
/* Allocate space */
yy = (int) (yymmdd / 10000);
mm = (yymmdd / 100) % 100;
dd = yymmdd % 100;
/* Is this a leap year? */
remainder = yy % 4;
leapyear = 0;
if (remainder == 0)
{
leapyear = 1;
}
/* Assign lengths to the months */
for (i=1; i<=12; i++)
{
month[i] = 31;
}
month[2] = leapyear ? 29 : 28;
month[4]=30;
month[6]=30;
month[9]=30;
month[11]=30;
/* Now just count up the days */
for (i=1; i<mm; i++)
{
jday += month[i];
}
jday += dd;
return (jday);
}
double airmass (double zeta)
{
/* cosine of zenith angle */
double costh;
double A;
costh = cos ((90-zeta) * 2 * M_PI / 360.0);
if (costh > 0)
{
A = 1.0/(costh + 0.50572 * exp (-1.6364 * log (6.07995 + zeta)));
return (A);
} /* if (costh... */
else
return (1e8);
} /* airmass (... */
#undef BARNARD_LANGLEY_C
/* end of file barnard_langley.c */
src/jim-plot.c 0 → 100644
View file @ 02287497
#include <stdio.h>
#include <math.h>
#include "jim-plot.h"
int jimmy_plot(plot_t *p)
{
int i, j;
FILE *f;
extern char *plot_fnm(void);
if ((f = fopen(plot_fnm(), "a")) == NULL) {
fprintf(stderr, "hey: can't open plot file\n");
return(-1);
}
fprintf(f, "%d\n", p->rows);
for (i=0; i<p->rows; i++) {
fprintf(f, "%d %.4f", i, p->am[i]);
for (j=0; j<p->cols; j++)
fprintf(f, " %.4f %d", log(p->points[j][i]), p->status_code[j][i]);
fputs("\n", f);
}
fclose(f);
return(1);
}
src/jim-plot.h 0 → 100644
View file @ 02287497
#ifndef JIM_PLOT_H
#define JIM_PLOT_H
typedef struct {
int cols, rows, type;
double *am, **points;
char *xlab, *ylab, *title,
**status_code;
} plot_t;
int jimmy_plot(plot_t *);
#endif
src/la.c 0 → 100644
View file @ 02287497
#define LA_C
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <errno.h>
#include "la.h"
#include "sunae.h"
#include "util.h"
#include "jim-plot.h"
#include "michalsky_langley.h"
/*
** 8/16/96 rs
**
** First try at making a 'standalone' langley analyzer.
** The intention here is to read in an ASCII file which contains
** several whitespace-separated records; each line of the file
** represents an observation, or record. Each field in each record
** should be a number representing various data relevant to solar
** energy measurement, particularly, the time and date of the observation,
** and the direct component of the sky at various wavelengths.
**
** The correspondence between fields in each record and what they
** actually represent is arbitrated by a configuration file. This
** config file is organized like a windows INI file...that is,
** each line looks like an assignment statement. The left side is
** the name of variable, the right side is the value of the variable,
** and in between is an equals ("=") sign. So,
YearCol = 1
** is a statement which indicates that the year of the observation
** is in column 1 of the input file (counting from 1).
** The job of the configuration file
** is to indicate to this program which fields are relevant, and what
** columns they occupy.
**
** The variables that must be assigned to are:
** Year (the year of the observation with or without century)
** Doy (the day of year, jan 1 == 1)
** Tod (the time of day of the observation)
** Latitude (Northern hemisphere positive)
** Longitude (West of Greenwich negative)
** HoursGMT (the number of hours to add to the given time to make it GMT)
** These last 3 aren't column numbers - just put the lat and lon
** as the arguments.
**
** The other important columns are those that contain direct normal
** radiation. This program can accommodate an arbitrary number of
** direct channels. The columns that these channels occupy are given with
** the "DirectCol" variable; for
** each Direct variable mentioned in the config file, a corresponding
** langley analysis is performed.
**
** This program reads the input file, and computes optical depth and
** extraterrestrial based on the data in it.
**
** -- The following has been changed...see 9/24/96 comment --
** It is TOTALLY up to the
** user to make sure that the data is partitioned into sensible
** units chronologically, and also that the data is calibrated,
** cosine corrected and so forth.
** -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
**
** The output consists of 2 numbers for each requested input channel:
** optical depth and extraterrestrial.
**
**
** 9/24/96 rs
** The input can now span an arbitrary time period. It will be split up
** into AM/PM groups and each tested separately.
**
**
** 2/11/97 rs
** Input data may be smoothed. This is provided because high sample rate
** data bounce around an awful lot and eliminate too many points.
** See smooth_input().
**
**
** 2/12/97 rs
** # Allow a bit of fudge when doing the first delta-difference test.
** # The la_cfg structure now has a new field called diff_slop which
** # can be set via command line or config file.
** Not really done yet! rs 2/19/97
** See 3/17/97 note below....
**
** 2/19/97 rs
** Implement "key wavelength" functionality.
**
**
** 2/27/97 rs
** "Users" want to have certain parameters in the filters be configurable
** at runtime. The first is the final lsfit standard deviation test; Lee
** had this set to 0.006, but make a config entry to allow for different
** values. Also add low and high airmass limits.
**
** 3/17/97 rs
** Officially put the 'cloud passage slop factor' in. The config file
** parameter name is "CloudSlop".
**
** 6/27/97 jim
** Add to to allow the user to enter dates in "julian" format rather
** than year,doy,tod.
**
** 6/27/97 rs
** Put a "-v" flag in that just prints the version number, which is
** now 1.0.
**
** 6/27/97 rs
** Also compute and print the number of records which are within the
** selected airmass range. v1.1
**
** 7/1/97 rs
** Make a special output style which is more parseable than human-readable.
** The command-line flag is "-p". v1.2
**
** 8/14/97 rs
** Work around the broken SunOS realloc. v 1.3
**
** 8/14/97 rs
** Try to speed it up on SunOS by reducing the number of calls to
** mktime() and sunae(). v 1.4
**
** 9/3/97 rs
** Install command-line flag and interface code so that we can call
** a plotting function for each langley interval. The actual plot
** function does not exist yet, but I can get ready for it. v 1.4.1
**
** 9/4/97 rs
** Fixed lots of problems; v1.4.2
**
** 1/29/98 jim
** Rob's gone but not forgotten. The use of mktime() in time_of() was
** misguided as mktime() returns a time_t corresponding to localtime.
** If your system runs on GMT (CUT, UTC) -- that is, your system time zone
** is set to 0 hours from GMT -- and your data was in GMT time,
** you may not have seen anything odd.
** On all other systems, however, you would have seen the langley report
** time given incorrectly. It would have been given as GMT + timezone
** which was really whacked. No effect on solar geometry or langley
** numbers though. v 1.4.3
**
** 3/5/98 jim
** Added automatic precision calculation to 'boring' output mode. This
** was needed in the case where the input irradiances were calibrated.
** This made the i0's quite small and the difference between the 1 AU
** normalized i0's and the 'local' i0's could not be seen.
** v 1.4.4
**
** 7/13/98 jim
** Something funky happening with getopt switch(). Under certain conditions
** a bad flag coupled with more args will screw up getopt's bookkeeping and
** generate a SEGV. Think I've made it a tad more bullet proof.
**
** 7/14/98 jim
** Found a file that would generate a log SING error. Added code in lsf that
** checks the incoming i0 for <= 0 so as to avoid SIGFPE messages.
** v 1.4.5
**
** 3/31/99 jim
** Added code to marry la with tu, the rsr unpacker. In fact, this code should
** make it fairly easy to marry la with any upstream data handler. I do this so
** that I don't need to upack to an ascii file and then suck the data back into
** la in order to do a langley.
**
** 8/22/2000 jim
** Late leap year fix in sunae.c
** v 1.4.6
*/
#define VERSION "1.4.6"
char *plot_fnm(void);
ae_pack *aeps;
time_t *tts;
/* tdh - moved hours_gmt, lat, lon to michalsky_langley.h */
void *my_realloc(void *, int);
/*
Note that if you're using la as a module and not a command line
version, you need to call the config_file() and langley() routines
from your calling code.
The calling code needs to include langley.h and must extern these
two globals:
rec_t *la_records = NULL;
int la_n_records = 0;
Memory needs to be allocated for la_records as well.
*/
#ifndef WITH_TU
int main(int argc, char **argv)
{
int args(int, char **),
config_file(void),
langley(void);
if (!args(argc, argv))
exit(1);
if (!config_file())
exit(1);
if (!langley())
exit(1);
return 0;
}
void do_version(char *exe)
{
fprintf(stderr, "%s version %s\n", exe, VERSION);
}
int args(int ac, char **av)
{
int c;
FILE *f;
extern int optind;
extern char *optarg;
void usage(char *),
do_version(char *);
arguments.cfg_fnm = NULL;
arguments.in_fnm = NULL;
arguments.debug = 0;
arguments.julian_out = 0;
arguments.version = 0;
arguments.smooth = 0;
arguments.parse_output = 0;
arguments.p_type = None;
while ((c = getopt(ac, av, "pvjhf:ds:x:")) != EOF) {
switch (c) {
case 'x':
if (strcmp(optarg, "ps") == 0)
arguments.p_type = PS;
else if (strcmp(optarg, "gif") == 0)
arguments.p_type = Gif;
else if (strcmp(optarg, "screen") == 0)
arguments.p_type = Screen;
else {
fprintf(stderr, "Error: -x: unknown plot type (%s)\n", optarg);
usage(av[0]);
return 0;
}
if ((f = fopen(plot_fnm(), "w")) == NULL) {
fprintf(stderr, "Error: Can't create plot file\n");
exit(1);
}
fclose(f);
break;
case 'j':
arguments.julian_out = 1;
break;
case 'p':
arguments.parse_output = 1;
break;
case 'v':
do_version(av[0]);
return 0;
case 'h':
usage(av[0]);
return 0;
case 's':
if (sscanf(optarg, "%d", &arguments.smooth) != 1) {
fprintf(stderr, "Error in '-s %s'; argument should "
"be an integer\n", optarg);
usage(av[0]);
return 0;
}
if (arguments.smooth <= 0 || arguments.smooth > 3600) {
fprintf(stderr, "Error in '-s %d'; argument should "
"be 1..3600\n", arguments.smooth);
usage(av[0]);
return 0;
}
break;
case 'f':
arguments.cfg_fnm = strdup(optarg);
break;
case 'd':
arguments.debug++;
break;
default:
return 0;
}
}
if (arguments.cfg_fnm == NULL)
arguments.cfg_fnm = strdup("la.cfg");
if (optind != ac - 1) {
fprintf(stderr, "No file to analyze! (optind = %d, ac = %d)\n", optind, ac);
return 0;
}
arguments.in_fnm = strdup(av[ac-1]);
return 1;
}
void usage(char *exe)
{
fprintf(stderr, "Usage: %s [-pdhv] [-x gr] [-s smooth] [-f cfgfile] infile\n", exe);
fputs("Where: '-p' produces more computer-readable output\n", stderr);
fputs(" '-j' used with -p, produces output similar to 'callang' w/julian dates\n", stderr);
fputs(" '-d' produces lots of debugging output\n", stderr);
fputs(" '-s secs' smooths the input to the "
"given # of seconds\n", stderr);
fputs(" '-f cfgfile' specifies the config file. Default is"
" la.cfg\n", stderr);
/*
fputs(" '-x gr' makes graphics output for device 'gr'\n", stderr);
fputs(" 'gr' can be 'ps', 'gif', or 'screen'\n", stderr);
*/
fputs(" '-v' prints version number\n", stderr);
fputs(" '-h' prints usage\n", stderr);
}
#endif /* WITH_TU */
int config_file(void)
{
char buf[512];
FILE *f;
int line,
chop_up(char *),
check_config(void);
if ((f = fopen(arguments.cfg_fnm, "r")) == NULL) {
fprintf(stderr, "Can't open config file '%s'\n", arguments.cfg_fnm);
return 0;
}
configuration.fields = NULL;
configuration.n_fields = 0;
configuration.julian = 0;
line = 0;
while (fgets(buf, 512, f) != NULL) {
line++;
if (!chop_up(buf)) {
fclose(f);
fprintf(stderr, "config file err: line %d\n", line);
return 0;
}
}
fclose(f);
return check_config();
}
int empty_line(char *s)
{
for ( ; *s; s++)
if (!isspace((int) *s))
return 0;
return 1;
}
int chop_up(char *s)
{
char *c;
int cleanup(field_t *);
if ((c = strchr(s, '#')) != NULL)
*c = '\0';
if ((c = strchr(s, '\n')) != NULL)
*c = '\0';
if (empty_line(s))
return 1;
if ((c = strchr(s, '=')) == NULL) {
fputs("config file: no '=' in line\n", stderr);
return 0;
}
*c = '\0';
if ((configuration.fields = my_realloc(configuration.fields,
(configuration.n_fields+1)*sizeof(field_t))) == NULL) {
fprintf(stderr, "alloc error in config file\n");
return 0;
}
configuration.fields[configuration.n_fields].checked = 0;
configuration.fields[configuration.n_fields].var = strdup(s);
configuration.fields[configuration.n_fields].val = strdup(c + 1);
if (!cleanup(configuration.fields + configuration.n_fields)) {
fputs("config file: empty field\n", stderr);
return 0;
}
configuration.n_fields++;
return 1;
}
int cleanup(field_t *f)
{
void trim_ws(char *);
if (empty_line(f->var))
return 0;
if (empty_line(f->val))
return 0;
trim_ws(f->var);
trim_ws(f->val);
return 1;
}
void trim_ws(char *s)
{
char *p;
/* remove trailing space */
for (p=strchr(s, '\0') - 1; p > s; p--)
if (isspace((int) *p))
*p = '\0';
else
break;
/* skip leading space */
for (p=s; isspace((int) *p); p++) {
}
/* copy byte-for-byte */
while (*s++ = *p++) {
}
}
int NonDirCols, n_columns = 0,
*xref = NULL;
/*
** go through 'configuration' and make sure that all required elements
** are present, and that they all make some kind of sense.
*/
int check_config(void)
{
int i, n_keys, temp,
required_double(char *, double *),
required_int(char *, int *),
direct_lookup(void),
int_lookup(char *, int *);
/* tdh - added the #ifndef */
#ifndef WITH_TU
if (!required_double("Lat", &lat))
return 0;
if (!required_double("Lon", &lon))
return 0;
if (!required_double("HoursGMT", &hours_gmt))
return 0;
#endif /* WITH_TU */
if (int_lookup("JulianCol", &temp)) {
NonDirCols = 1;
if ((xref = malloc(sizeof(double))) == NULL) {
fprintf(stderr, "Alloc fail: 1 double in check_config\n");
exit(1);
}
*xref = temp;
configuration.julian = 1;
if(int_lookup("YearCol", &temp) ||
int_lookup("DoyCol", &temp) ||
int_lookup("TodCol", &temp)) {
fprintf(stderr, "JulianCol cannot coexist with other YearCol, DoyCol or TodCol.\n");
exit(1);
}
}
else {
NonDirCols = 3;
if ((xref = malloc(3*sizeof(double))) == NULL) {
fprintf(stderr, "Alloc fail: 3 doubles in check_config\n");
exit(1);
}
if (!required_int("YearCol", xref))
return 0;
if (!required_int("DoyCol", xref + 1))
return 0;
if (!required_int("TodCol", xref + 2))
return 0;
}
n_columns = NonDirCols;
n_keys = 0;
while ((i = direct_lookup()) > 0) {
if ((xref = my_realloc(xref, (n_columns+1)*sizeof(int))) == NULL) {
fprintf(stderr, "too many directs-alloc fail at %d\n",
n_columns+1);
exit(1);
}
xref[n_columns] = i;
n_columns++;
if (i > 1000)
n_keys++;
}
if (n_columns == NonDirCols) {
fprintf(stderr, "Config: There are no direct columns!\n");
free(xref);
return 0;
}
if (n_keys > 1) {
fprintf(stderr, "Config: A maximum of 1 key fields is allowed\n");
free(xref);
return 0;
}
/* config file column numbers are 1-based; change to 0 now. */
for (i=0; i<n_columns; i++)
xref[i]--;
if (arguments.debug > 0) {
for (i=0; i<configuration.n_fields; i++)
printf("'%s' '%s'\n", configuration.fields[i].var,
configuration.fields[i].val);
printf("Lat %.3f\n", lat);
printf("Lon %.3f\n", lon);
switch (NonDirCols) {
case 3:
printf("Year %d\n", xref[0]);
printf("Doy %d\n", xref[1]);
printf("Tod %d\n", xref[2]);
break;
case 1:
printf("Julian %d\n", xref[0]);
break;
default:
fprintf(stderr, "Unknown number of NonDirCols: %d\n", NonDirCols);
exit(1);
}
printf("Directs at");
for (i=NonDirCols; i<n_columns; i++)
if (xref[i] >= 1000)
printf(" %d Key", xref[i] - 1000);
else
printf(" %d", xref[i]);
puts("");
}
return 1;
}
int int_lookup(char *key, int *x)
{
char *s, *lookup(char *);
if ((s = lookup(key)) == NULL)
return 0;
if (sscanf(s, "%d", x) != 1) {
fprintf(stderr, "Config file: Failed to convert value of '%s' to"
"an integer\n", key);
return 0;
}
return 1;
}
int double_lookup(char *key, double *x)
{
char *s, *lookup(char *);
if ((s = lookup(key)) == NULL)
return 0;
if (sscanf(s, "%lf", x) != 1) {
fprintf(stderr, "Config file: Failed to convert value of '%s' to"
"a double\n", key);
return 0;
}
return 1;
}
int direct_lookup(void)
{
int i, n,
ig_strcmp(char *, char *);
char val[1024], *sp[64];
for (i=0; i<configuration.n_fields; i++)
if (strcmp(configuration.fields[i].var, "DirectCol") == 0 &&
!configuration.fields[i].checked) {
configuration.fields[i].checked++;
strncpy(val, configuration.fields[i].val, 1023);
n = split(val, sp, 64, "");
if (sscanf(sp[0], "%d", &i) != 1) {
fprintf(stderr, "config file: direct column err: "
"what is '%s'?\n", sp[0]);
exit(1);
}
if (i <= 0) {
fprintf(stderr, "config file: direct column bad: "
"'%d'?\n", i);
exit(1);
}
if (n > 0 && ig_strcmp(sp[1], "key") == 0)
i += 1000; /* cheesey but true */
return i;
}
return -1;
}
char *lookup(char *key)
{
int i;
for (i=0; i<configuration.n_fields; i++)
if (strcmp(key, configuration.fields[i].var) == 0)
return configuration.fields[i].val;
/* fprintf(stderr, "Failed to find key '%s' in config file", key); */
return NULL;
}
/*
** read the input file and convert all the required columns
*/
int langley(void)
{
int i, input_loop(void);
void smooth_input(int),
times_sunaes(),
langley_analyze(void);
#ifndef WITH_TU
if (!input_loop())
return 0;
#endif
times_sunaes();
smooth_input(arguments.smooth);
langley_analyze();
if (la_n_records != 0) {
for (i=0; i<la_n_records; i++)
free(la_records[i].direct);
free(la_records);
}
free(aeps);
free(tts);
return 1;
}
int input_loop(void)
{
FILE *f;
char buf[512], *sp_buf[512], *p;
int i, lines, julian_2_ydm();
rec_t *rec;
extern int errno;
static int cur_size = 0;
if ((f = fopen(arguments.in_fnm, "r")) == NULL)
{
fprintf(stderr, "Can't open input file '%s'\n", arguments.in_fnm);
return 0;
}
lines = 0;
while (fgets(buf, 512, f) != NULL) {
lines++;
if ((p = strchr(buf, '\n')) != NULL)
*p = '\0';
if (split(buf, sp_buf, 512, "") < n_columns) {
fprintf(stderr, "Input file: line %d: too few fields\n", lines);
lines = 0;
break;
}
if (la_n_records == cur_size) {
if (cur_size == 0)
cur_size = 64;
else
cur_size *= 2;
if ((la_records=my_realloc(la_records,cur_size*(sizeof(rec_t)))) == NULL) {
fprintf(stderr, "input: alloc failure at line %d\n", lines);
exit(1);
}
}
rec = la_records + la_n_records;
errno = 0;
if(configuration.julian) {
if(!julian_2_ydm(rec, strtod(sp_buf[xref[0]], NULL)))
exit(0);
}
else {
rec->year = atoi(sp_buf[xref[0]]);
rec->doy = atoi(sp_buf[xref[1]]);
rec->tod = strtod(sp_buf[xref[2]], NULL);
}
rec->n_direct = n_columns - NonDirCols;
if ((rec->direct = malloc(rec->n_direct*sizeof(double))) == NULL) {
fprintf(stderr, "alloc fail: directs #%d\n", la_n_records);
exit(1);
}
for (i=NonDirCols; i<n_columns; i++)
rec->direct[i - NonDirCols] =
strtod(sp_buf[xref[i] % 1000], NULL);
if (errno != 0) {
fprintf(stderr, "input: conversion error line %d\n", lines);
lines = 0;
break;
}
la_n_records++;
}
fclose(f);
return lines;
}
void times_sunaes(void)
{
int i;
obs_time_t o;
if ((aeps = malloc(la_n_records*sizeof(ae_pack))) == NULL) {
fprintf(stderr, "times_sunaes: malloc fail (%d) ae_pack\n", la_n_records);
exit(1);
}
if ((tts = malloc(la_n_records*sizeof(time_t))) == NULL) {
fprintf(stderr, "times_sunaes: malloc fail (%d) time_t\n", la_n_records);
exit(1);
}
for (i=0; i<la_n_records; i++) {
aeps[i].is_tst = 0;
aeps[i].lat = lat;
aeps[i].lon = lon;
aeps[i].year = la_records[i].year;
aeps[i].doy = la_records[i].doy;
aeps[i].hour = la_records[i].tod + hours_gmt;
if (aeps[i].hour > 24.0) {
aeps[i].hour -= 24.0;
if (++aeps[i].doy > 366) {
aeps[i].doy = 1;
aeps[i].year++;
}
}
if (aeps[i].hour < 0.0) {
aeps[i].hour += 24.0;
if (--aeps[i].doy < 1) {
aeps[i].doy = 365;
aeps[i].year--;
}
}
(void) sunae(aeps + i);
o.year = aeps[i].year;
o.doy = aeps[i].doy;
o.gmt = aeps[i].hour;
tts[i] = time_of(&o);
}
}
void langley_analyze(void)
{
int i, j, k, l, n, p;
void print_records(int, int, int),
free_lrv(la_return *, int),
emit_lang(la_return *l, la_cfg *cfg, la_obs *);
double dt, partition(int, int *, int *);
la_cfg la_config;
la_obs *obs;
la_return *la_ret;
la_config.lat = lat;
la_config.lon = lon;
/* - tdh if (!double_lookup("LsfitSD", &la_config.lsfit_sd)) */
la_config.lsfit_sd = gLSFITSD;
/* - tdh if (!double_lookup("LowAM", &la_config.low_am)) */
la_config.low_am = LOW_AM;
/* - tdh if (!double_lookup("HighAM", &la_config.high_am)) */
la_config.high_am = HIGH_AM;
/* - tdh if (!double_lookup("OutLimit", &la_config.out_limit)) */
la_config.out_limit = OUTLIMIT;
/* - tdh if (!double_lookup("FracPts", &la_config.frac_pts)) */
la_config.frac_pts = FRACPTS;
/* - tdh if (!double_lookup("CloudSlop", &la_config.diff_slop)) */
la_config.diff_slop = CLOUDSLOP;
/* - tdh n = n_columns - NonDirCols; */
n = NCHANNELS;
la_config.n_data = n;
if ((la_config.ana = malloc(n*sizeof(ana_t))) == NULL) {
fprintf(stderr, "can't alloc %d ana_t\n", n);
exit(1);
}
for (i=0; i<n; i++)
/* - tdh if (xref[i + NonDirCols] >= 1000) */
if (i == KEY_CHN)
la_config.ana[i] = Key_Ana;
else
la_config.ana[i] = Do_Ana;
for (p=0; (dt = partition(p, &j, &k)) > 0.0; p++) {
n = k - j + 1;
// When I correct for azimuth, we often end up with orphaned
// partitions, so demand at least three points to use this
// partition
if (n < 3) { continue;}
if ((obs = malloc(n*sizeof(la_obs))) == NULL) {
fprintf(stderr, "can't alloc %d la_obs\n", la_n_records);
exit(1);
}
for (i=j; i<=k; i++) {
l = i - j;
obs[l].obs_time.year = aeps[i].year;
obs[l].obs_time.doy = aeps[i].doy;
obs[l].obs_time.gmt = aeps[i].hour;
obs[l].obs_time.tt = tts[i];
obs[l].data = la_records[i].direct;
//printf ("la_records ... i=%d\n");
//printf (" direct=%lf\n", *la_records[i].direct);
}
if (arguments.debug > 0)
print_records(p, j, k);
la_config.n_obs = n;
la_ret = lang_ana(&la_config, obs);
/* - tdh This is not needed by the dmf/
emit_lang(la_ret, &la_config, obs);
*/
free_lrv(la_ret, la_config.n_data);
free(obs);
}
free(la_config.ana);
}
/*
** A "partition" of the input data is a group of records which are all
** either before or after local noon for a single day. This function
** returns the indeces of the first and last record which correspond
** to partition number 'part,' counting from 0. If there is no partition
** with that number, return 0; otherwise return 1.
*/
double partition(int part, int *first, int *last)
{
static int done = 0, *p;
double rt;
int i, *make_partition_list(void);
if (!done) {
p = make_partition_list();
done = 1;
}
for (i=0; i<la_n_records; i++)
if (p[i] == part) {
*first = i;
break;
}
if (i == la_n_records)
return 0.0;
while (i < la_n_records && p[i] == part)
*last = i++;
rt = 1000.0*la_records[*first].year + la_records[*first].doy;
rt += la_records[*first].tod/24.0;
return rt;
}
int *make_partition_list(void)
{
int i, *p, cp;
int diff_half_day(int, int);
if ((p = malloc(la_n_records*sizeof(int))) == NULL) {
fprintf(stderr, "make_partition_list: can't allocate vector\n");
exit(1);
}
cp = 0;
p[0] = 0;
for (i=1; i<la_n_records; i++) {
if (diff_half_day(i-1, i))
cp++;
p[i] = cp;
}
return p;
}
/*
** if t1 and t2 are not in the same half of the same day, return 1.
** Otherwise, return 0.
*/
int diff_half_day(int t1, int t2)
{
time_t dt;
dt = tts[t2] - tts[t1];
if (dt < 0)
dt = -dt;
if (dt > 43200)
return 1;
// TRS: Believe it or not, this works to find noon and midnight
// boundaries for both the northern and southern hemispheres. In the
// north, az increases throughout the day, from 360/0 at midnight to
// 180 at noon. In the southern hemisphere, az decreases throughout
// the day, from 180 at midnight to 0/360 at noon. In either case, and
// for both half-day boundaries, this returns true if and only iff t1
// and t2 are on opposite sides of the half-day.
if (aeps[t1].az < 180.0)
return aeps[t2].az >= 180.0;
else
return aeps[t2].az < 180.0;
}
void print_records(int part, int first, int last)
{
int i, j;
if (arguments.debug > 0) {
printf("---- Partition %d has %d records ----", part, last-first+1);
if (arguments.debug > 1) {
puts(":");
for (i=first; i<=last; i++) {
printf("%d %d %.2f", la_records[i].year, la_records[i].doy,
la_records[i].tod);
for (j=0; j<la_records[i].n_direct; j++)
printf(" %.1f", la_records[i].direct[j]);
puts("");
}
} else
puts(".");
}
}
/*
** make a pass through the input data and make each data point the
** average of the nearest several points surrounding it. the size
** of the "window" is determined by the 'size' parameter.
*/
void smooth_input(int size)
{
int i, j;
double *list;
void smooth(time_t *, double *, int, int);
if (size < 1)
return;
if (la_n_records < 1)
return;
if ((list = malloc(la_n_records*sizeof(double))) == NULL) {
fprintf(stderr, "malloc fail: smoothing input (%d doubles)\n",
la_n_records);
exit(1);
}
for (j=0; j<la_records[0].n_direct; j++) {
for (i=0; i<la_n_records; i++)
list[i] = la_records[i].direct[j];
smooth(tts, list, size, la_n_records);
for (i=0; i<la_n_records; i++)
la_records[i].direct[j] = list[i];
}
free(list);
}
/*
** The input consists of 2 lists of length 'n'; 'tt' is the time of
** each observation, and 'pts' are the data values. 'size' is the half size
** of the window around which to smooth each point.
*/
void smooth(time_t *tt, double *pts, int size, int n)
{
int i, j, first, count, smoothed;
time_t t1, t2;
double sum;
smoothed = 0;
first = 0;
for (i=0; i<n; i++) {
t1 = tt[i] - size;
t2 = tt[i] + size;
sum = 0.0;
count = 0;
for (j=first; j<n; j++)
if (tt[j] < t1)
continue;
else if (tt[j] > t2)
break;
else {
sum += pts[j];
count++;
if (count == 1)
first = j;
}
if (count > 0) /* should always be... */
pts[i] = sum/count;
if (count > 1)
smoothed++;
}
if(arguments.debug)
fprintf(stderr, "%d out of %d points were smoothed\n", smoothed, n);
}
int ig_strcmp(char *s1, char *s2)
{
if (s1 == NULL)
if (s2 == NULL)
return 0;
else
return 1;
else
if (s2 == NULL)
return 1;
while (*s1 != '\0') {
if (toupper(*s1) != toupper(*s2))
return 1;
s1++;
s2++;
}
return *s2 != '\0';
}
void no_cfg_entry(char *s)
{
fprintf(stderr, "Config: Required config entry '%s' not found\n", s);
}
int required_double(char *key, double *v)
{
int double_lookup(char *, double *);
if (!double_lookup(key, v)) {
no_cfg_entry(key);
return 0;
} else
return 1;
}
int required_int(char *key, int *v)
{
int int_lookup(char *, int *);
if (!int_lookup(key, v)) {
no_cfg_entry(key);
return 0;
} else
return 1;
}
#define DAYS_1900_1970 25568L /* days from 1/1/00 to 1/1/1970 */
int julian_2_ydm(rec_t *rec, double jdays) {
long lround(double);
struct tm *ts;
time_t secs;
jdays -= (double) DAYS_1900_1970;
if(jdays < 0.0) {
fprintf(stderr, "Bad julian date: %.5f (must be >= %ld)\n", jdays, DAYS_1900_1970);
return 0;
}
secs = (time_t) lround( (jdays * 3600.0 * 24.0) );
ts = gmtime(&secs);
rec->year = ts->tm_year;
rec->doy = ts->tm_yday + 1;
rec->tod = (double) ts->tm_hour + ts->tm_min/60.0 + ts->tm_sec/3600.0;
return 1;
}
double unix_2_julian(double unix_time) {
double jdays = (double) DAYS_1900_1970;
jdays += unix_time / 86400.0;
return(jdays);
}
/* prolly won't work for negative numbers */
long lround(double d)
{
long retval;
double remainder;
retval = (long) d;
remainder = d - (double) retval;
return((remainder >= 0.5) ? retval+1 : retval);
}
void print_dates(time_t t1, time_t t2)
{
char buf[512];
struct tm *stm;
double unix_2_julian();
if (!arguments.parse_output) {
fputs("Time span: ", stdout);
stm = gmtime(&t1);
strftime(buf, 512, "%m/%d/%y %H:%M:%S to ", stm);
fputs(buf, stdout);
stm = gmtime(&t2);
strftime(buf, 512, "%m/%d/%y %H:%M:%S", stm);
fputs(buf, stdout);
fputs("\n", stdout);
}
else if(arguments.julian_out) {
/* if the -j flag was given, print out the midpoint of the
airmass period in julian time and don't stick a CR in there */
printf("%.5f ", unix_2_julian((double)t1 + (double) (t2-t1)/2.0));
}
else {
stm = gmtime(&t1);
strftime(buf, 512, "%m/%d/%y %H:%M:%S ", stm);
fputs(buf, stdout);
stm = gmtime(&t2);
strftime(buf, 512, "%m/%d/%y %H:%M:%S", stm);
fputs(buf, stdout);
fputs("\n", stdout);
}
}
void emit_lang(la_return *l, la_cfg *cfg, la_obs *obs)
{
void boring_output(la_return *, la_cfg *, la_obs *),
pretty_output(la_return *, la_cfg *, la_obs *),
plot_output(la_return *, la_cfg *, la_obs *);
if (arguments.parse_output)
boring_output(l, cfg, obs);
else
pretty_output(l, cfg, obs);
if (arguments.p_type != None)
plot_output(l, cfg, obs);
}
void boring_output(la_return *l, la_cfg *cfg, la_obs *obs)
{
int i, n;
void print_dates(time_t, time_t);
char ioprec[10], format_str[256], *fmt="%.5f %.5f %d/%d/%d\n";
/* if i0 is calibrated, then it's going to be smaller and we
need more precision
*/
if(l[0].i0 < 10)
strcpy(ioprec, "%.4f");
else if(l[0].i0 < 100)
strcpy(ioprec, "%.3f");
else
strcpy(ioprec, "%.2f");
sprintf(format_str, "%s %s %s", ioprec, ioprec, fmt);
for (n=0, i=0; n == 0 && i<cfg->n_data; i++)
if (l[i].op_depth != 0.0 && l[i].i0 != 0.0 && l[i].n != 0)
n++;
if (n == 0)
return;
if(!arguments.julian_out)
print_dates(l->first_am, l->last_am);
for (i=0; i<cfg->n_data; i++)
if (l[i].op_depth != 0.0 && l[i].i0 != 0.0 && l[i].n != 0) {
if(arguments.julian_out)
print_dates(l->first_am, l->last_am);
printf("%d ", i+1);
printf(format_str,
l[i].i0, l[i].i0_sdist, -l[i].op_depth,
l[i].sd, l[i].n, l->n_am, cfg->n_obs);
}
}
void pretty_output(la_return *l, la_cfg *cfg, la_obs *obs)
{
int i;
void print_dates(time_t, time_t);
print_dates(obs[0].obs_time.tt, obs[cfg->n_obs-1].obs_time.tt);
for (i=0; i<cfg->n_data; i++) {
printf("Channel %3d: ", i+1);
if (l[i].op_depth == 0.0 && l[i].i0 == 0.0)
printf("%7.2f %7.2f %9.5f %8.5f %5d/%d/%d No points to regress\n",
0.0, 0.0, 0.0, 0.0, 0, l->n_am, cfg->n_obs);
else if (l[i].n == 0)
printf("%7.2f %7.2f %9.5f %8.5f %5d/%d/%d All points filtered.\n",
l[i].i0, l[i].i0_sdist, -l[i].op_depth,
l[i].sd, l[i].n, l->n_am, cfg->n_obs);
else
printf("%7.2f %7.2f %9.5f %8.5f %5d/%d/%d\n",
l[i].i0, l[i].i0_sdist, -l[i].op_depth,
l[i].sd, l[i].n, l->n_am, cfg->n_obs);
}
}
void *my_realloc(void *p, int n)
{
if (p == NULL)
return malloc(n);
else
return realloc(p, n);
}
void plot_output(la_return *l, la_cfg *cfg, la_obs *obs)
{
plot_t p;
int i, j;
int status;
if (l[0].n_am == 0)
return;
p.type = arguments.p_type;
p.rows = l[0].n_am;
p.cols = cfg->n_data;
p.xlab = strdup("Airmass");
p.ylab = strdup("LOG(Direct)");
p.title = strdup("Langley Analysis");
if ((p.am = malloc(p.rows*sizeof(double))) == NULL) {
fprintf(stderr, "plot_output: malloc (1) fail\n");
exit(1);
}
if ((p.points = malloc(p.cols*sizeof(double *))) == NULL) {
fprintf(stderr, "plot_output: malloc (2) fail\n");
exit(1);
}
if ((p.status_code = malloc(p.cols*sizeof(char *))) == NULL) {
fprintf(stderr, "plot_output: malloc (3) fail\n");
exit(1);
}
for (i=0; i<p.cols; i++) {
if ((p.points[i] = malloc(p.rows*sizeof(double))) == NULL) {
fprintf(stderr, "plot_output: malloc (4) fail\n");
exit(1);
}
if ((p.status_code[i] = malloc(p.rows)) == NULL) {
fprintf(stderr, "plot_output: malloc (5) fail\n");
exit(1);
}
}
for (i=0; i<l[0].n_am; i++) {
p.am[i] = l[0].x[i];
for (j=0; j<cfg->n_data; j++) {
p.points[j][i] = l[j].y[i].y;
p.status_code[j][i] = l[j].y[i].flag;
}
}
status = jimmy_plot(&p);
if (status < 0)
{
return;
}
for (i=0; i<p.cols; i++) {
free(p.points[i]);
free(p.status_code[i]);
}
free(p.am);
free(p.points);
free(p.status_code);
}
void free_lrv(la_return *lrv, int n)
{
int i;
void free_lrv_1(la_return *);
if (lrv == NULL)
return;
for (i=0; i<n; i++)
free_lrv_1(lrv + i);
free(lrv);
}
void free_lrv_1(la_return *l)
{
if (l->n_am != 0) {
free(l->x);
free(l->y);
}
}
char *plot_fnm(void)
{
switch (arguments.p_type) {
case None:
return NULL;
case PS:
return "la-out.ps";
case Gif:
return "la-out.gif";
case Screen:
return "la-out.screen";
}
}
#undef LA_C
src/la.h 0 → 100644
View file @ 02287497
#ifndef LA_H
#define LA_H
#include <time.h>
#include "michalsky_langley.h"
/* tdh moved ana_t to michalsky_langley.h */
/* tdh moved la_cfg to michalsky_langley.h */
/* tdh moved obs_time_t to michalsky_langley.h */
/* tdh moved la_obs to michalsky_langley.h */
/* tdh - moved pt_t to michalsky_langley.h */
/* tdh - moved la_return to michalsky_langley.h */
/*
** Do langley analysis for a set of observations. The input records should
** be partitioned into whatever conditions are required.
** 'cfg' points to a description of the input data
** 'obs' points to a list of cfg->n_obs input records
**
** the return value is a pointer to a list of optical depth/i0 pairs; the
** number of pairs in the list is determined by how many 'columns' have
** langley analysis requested in cfg->ana.
*/
/* tdh - moved lang_ana(..) to michalsky_langley.h */
time_t time_of(obs_time_t *);
#ifndef M_DTOR
#define M_DTOR 0.0174532925199433
#define M_RTOD 57.2957795130823230
#define M_2PI 6.2831853071795862320E0
#define M_HTOR 0.2617993877991494
#define M_RTOH 3.8197186342054881
#define M_HTOD 15.0
#define M_DTOH 0.0666666666666667
#endif
#endif
src/langley.c 0 → 100644
View file @ 02287497
/*******************************************************************************
* COPYRIGHT (C) 2001 Battelle Memorial Institute.
* All Rights Reserved. (Now, with this out of the way...)
*
* RCS INFORMATION:
* $RCSfile: langley.c,v $
* $Revision: 1.40 $
* $Author: koontz $
* $Locker: $
* $Date: 2012/02/06 19:24:24 $
* $State: Exp $
* $Name: $
* $Id: langley.c,v 1.40 2012/02/06 19:24:24 koontz Exp $
*
* COPYRIGHT (C) 2001 Battelle Memorial Institute. All Rights Reserved
*
* FUNCTIONS IN THIS FILE:
* DATA* ComputeNewData (DATA *D)
*
* DESIGN:
* <Discuss overall design/purpose of class.>
*
* $Log: langley.c,v $
* Revision 1.40 2012/02/06 19:24:24 koontz
* *** empty log message ***
*
* Revision 1.39 2012/01/10 17:41:37 koontz
* added debug (temporary) to help figure out a problem in production.
*
* Revision 1.38 2012/01/04 19:36:35 koontz
* Added some debug print so we can figure out why E31 won't run in production.
*
* Revision 1.37 2011/12/29 18:56:50 koontz
* *** empty log message ***
*
* Revision 1.36 2010/02/09 14:05:52 koontz
* fixed a problem with mfrsr/nimfr
*
* Revision 1.35 2010/01/22 20:45:11 koontz
* Uncalibrate langley results so that frequent calibration of MFRSR or NIMFR will not adversely affect results.
*
* Revision 1.34 2009/09/30 15:25:47 koontz
* Dropping computation of the angstrom exponent field, so we don't need atmospheric pressure as an input.
*
* Revision 1.33 2009/02/03 20:52:11 koontz
* *** empty log message ***
*
* Revision 1.32 2008/09/25 15:41:12 koontz
* fix "pfound" when altitude is used.
*
* Revision 1.31 2008/08/22 19:10:14 koontz
* change how atmos pres is stored.
*
* Revision 1.30 2008/08/20 17:14:36 koontz
* fix to getPresFraction so that if it looks harder for pfrac
*
* Revision 1.29 2008/08/13 20:48:38 koontz
* add some logging info on what happens in getPresFraction so that we can (maybe) figure out what is happening in production.
*
* Revision 1.28 2008/07/21 19:33:04 koontz
* modified logic so that one can specify a different value for lsfitsd. Prior to this release, the default
* for this value was 0.006. But, runs for difficult sites (Niamey, for example) show a value of 0.018 produces
* better results. For SGP, this value of 0.018 seems to produce better results too. We added a new parameter to
* the commandline, so if one wants a different value, specify -x x.xxx, and that x.xxx value will be used.
*
* And, thanks to Robin Perez' keen eyes, we discovered a bug in how the surface pressure is calculated. This release fixes that bug.
*
* Revision 1.27 2008/01/16 20:38:22 koontz
* *** empty log message ***
*
* Revision 1.26 2008/01/16 20:23:01 koontz
* mods to fix NaN in getPresFraction, and to fix NaN/Inf in michalsky angstrom exponent.
*
* Revision 1.25 2007/11/19 20:01:19 koontz
* *** empty log message ***
*
* Revision 1.24 2007/11/19 19:33:17 koontz
* *** empty log message ***
*
* Revision 1.23 2007/11/15 13:41:46 koontz
* fixes to eliminate NaN/Inf/uninitialized values.
*
* Revision 1.22 2007/02/08 14:17:04 koontz
* *** empty log message ***
*
* Revision 1.21 2006/12/21 19:00:24 koontz
* *** empty log message ***
*
* Revision 1.20 2006/12/06 14:26:30 koontz
* cleaned up some non-used options, made some compiler warnings go away, etc.
*
* Revision 1.19 2006/12/05 20:36:58 koontz
* changes to implement input datastreams to header, misc. changes to find
* channel wavelengths as floats, not integers, changes to pressure fraction calculations
* (detect multiple units from observations and calculate properly.)
*
* Revision 1.18 2006/11/16 13:24:04 koontz
* *** empty log message ***
*
* Revision 1.17 2005/08/19 18:38:05 koontz
* Added versioning info to the output netCDF files, per BCR 1073.
* Some tidying up, and removal of A0 and A1 MFR input.
*
* Revision 1.16 2005/02/15 14:51:57 koontz
* Still another attempt to get quicklooks to work in production.
*
* Revision 1.15 2005/02/02 15:03:36 koontz
* adding info on command-line args so we can (maybe) figure
* out why quicklooks aren't being done in production.
*
* Revision 1.14 2004/11/08 16:09:56 halter
* changed presfraction to presFraction
*
* Revision 1.13 2002/10/22 17:12:41 halter
* fixed wavelength attr.
*
* Revision 1.12 2002/07/22 16:12:50 halter
* Getting ready for release.
*
* Revision 1.11 2002/04/18 16:20:46 halter
* Save before changes to make it work at twp.
*
* Revision 1.10 2001/11/29 18:08:42 halter
* barnard and michalsky now in same file.
*
* Revision 1.9 2001/11/08 17:57:59 halter
* Check in before merging the outputs.
*
* Revision 1.8 2001/09/10 17:41:34 halter
* *** empty log message ***
*
* Revision 1.7 2001/08/02 17:55:21 halter
* a good time to check in the code.
*
* Revision 1.6 2001/07/11 16:39:29 halter
* *** empty log message ***
*
* Revision 1.5 1998/12/09 00:24:12 shippert
* Ported to Y2K and BW3.3/Zebra_ASK. Also use bw_abort() instead of the
* homebrewed abort_vap(), and changed our matherr to a matherr_hook().
* Finally, added a hook to change platform name based on -S or -D cmdline
* options, thereby eliminating the annoying links to Slangley and Dlangley.
*
* Revision 1.4 1998/02/19 17:42:38 d3a230
* Added a check to see if our last langley plot was complete, i.e. if the
* last sample fell outside the airmass range. If not, then we send a
* warning, and we may wish to rerun when new data is available. Also added a
* matherr function, so that matherrs will output to the Eventlogger.
* (Probably will have to remove this in the new version of BW).
*
* Revision 1.3 1997/07/23 18:43:00 d3a230
* Played around with the breakout condition for our regression loop. Now we
* breakout if we don't reject any more points, which doesn't change anything
* but is more aesthetically pleasing than going the whole 20 runs.
*
* Revision 1.2 1997/07/22 16:18:46 d3a230
* Added Angstrom exponents and fixed some minor bugs. This version should be
* ready for release.
*
* Revision 1.1 1996/12/06 17:00:47 d3a230
* Initial revision
*
******************************************************************************/
#define LANGLEY_C
/****** System Includes ******/
#include <stdio.h>
#include <string.h>
#include <math.h>
/* #include <varargs.h> */
#include <stdarg.h>
#include <stdlib.h>
#include <sys/types.h>
#include <unistd.h>
#include <sys/stat.h>
#include <dirent.h>
#include <ctype.h>
#include "netcdf.h"
/****** Zebra Includes ******/
/* #include "config.h"
* #include "defs.h"
* #include "message.h"
* #include "timer.h"
* #include "DataStore.h"
*/
/****** BW Includes ******/
#define CUSTOM_HOOKS
/* #include "bw_main.h"
* #include "bw_retriever.h"
* #include "bw_unPacker.h"
* #include "bw_maker.h"
*/
#include "langley_retriever.h"
/****** RDSS User Interface ******/
/*#include <ui_error.h>*/
/****** Application Includes ******/
#include "langley.h"
/* Required for RCS */
static char *rcsid="$Id: langley.c,v 1.40 2012/02/06 19:24:24 koontz Exp $";
static char *rcsstate="$State: Exp $";
static char *rcsdate="$Date: 2012/02/06 19:24:24 $";
char gPlatNames[NPRESPLATS][64];
void custom_define_hooks(void);
/* global data used for command line arguments and hooks to modify input data */
char platform[20];
char datalevel[20];
#define FALSE 0
#define TRUE 1
/* globals used for quicklooks */
int rundate = -9999;
int plot_only_flag = FALSE;
int use_x_device = FALSE;
int quicklook_flag = TRUE;
int save_flag = FALSE;
int dev_flag = FALSE;
long startime;
long etime;
char gSite[20];
char gFacility[20];
void getWavelengthAttr ();
DATA* ComputeNewData (DATA *D)
/*****************************************************************************
* Input parameters:
*
* Returned value:
*
* Called from: main
*
* Functions called:
*
* Description:
* This function is invoked by BW for every application mode (ingest,
* measurement, and delivery). For data pass-through, nothing needs to
* be added. Simply return the input parameter ("D" in this template).
*
*****************************************************************************/
{
DATA *newD;
/* function declarations */
int chechQCFlag (DATA *);
int getStartEndTime (DATA *, long *, long *);
int getStartEndTime3 (DATA *, long *, long *);
DATA *barnard_langley (DATA *, float *);
int michalsky_langley (DATA *, DATA *, float *);
int writeWavelengthAttr (DATA *);
char *head_id;
char *logger_id;
int o;
int status;
void checkQCFlag();
float gNomCal[7];
/* First thing, store global variables gSite and gFacility */
sprintf (gSite, "%s", dsproc_get_site());
sprintf (gFacility, "%s", dsproc_get_facility());
strcpy (&platform[0], dsproc_get_name());
/* determine which platform to use */
/* We did that in parse_commandline_options ... so, we will always want to use "p = 0" */
B1 = 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));
}
/* assume all obs are at the same location */
/* set global position vars. */
lat = (double) D->subLoc[B1][0][0].l_lat;
lon = (double) D->subLoc[B1][0][0].l_lon;
alt = (double) D->subLoc[B1][0][0].l_alt;
/* aaa */
/* Get head id and unit ID from b1 data */
head_id = (char *) bw_GetAttrValue (D,
B1,
0,
-1,
"head_id",
NULL,
NULL);
logger_id = (char *) bw_GetAttrValue (D,
B1,
0,
-1,
"logger_id",
NULL,
NULL);
getWavelengthAttr (D);
/* bbb */
/* printf ("lat=%lf ... lon=%lf ... alt=%lf\n", lat, lon, alt); */
checkQCFlag (D);
/* We will need to get the "rundate". This used to
* come in via the "-d" command line argument.
* ADI uses the -b concept, so I've deactivated
* the parse_commandline()
*
* Try using the retriever_GetStartTime() function instead
*/
if ((status = getStartEndTime (D, &startime, &etime)) < 0) {
bw_exit(NULL, "Problem getting start or end time; continuing to next processing interval");
}
/* Get head id and unit ID from b1 data */
head_id = (char *) bw_GetAttrValue (D,
B1,
0,
-1,
"head_id",
NULL,
NULL);
logger_id = (char *) bw_GetAttrValue (D,
B1,
0,
-1,
"logger_id",
NULL,
NULL);
getWavelengthAttr (D);
if ((newD = barnard_langley (D, (float *) &gNomCal)) == NULL) {
bw_exit(newD, "Problem in Barnard; continuing to next processing interval");
}
if ( (status = michalsky_langley (D, newD, (float *) &gNomCal)) < 0) {
bw_exit(newD, "Problem in Michalsky; continuing to next processing interval");
};
writeWavelengthAttr (newD);
for (o=0; o < newD->nObs[LANGLEY]; o++)
{
bw_SetAttrValue(newD,
LANGLEY,
o,
-1,
DCT_String,
1,
"head_id",
(void *) head_id);
bw_SetAttrValue(newD,
LANGLEY,
o,
-1,
DCT_String,
1,
"logger_id",
(void *) logger_id);
}
/* Add global metadata describing input data
* per BCR 1073
*/
return (newD);
} /* ComputeNewData */
char * custom_GetVersion()
/****************************************************************************
* This function is called by dc_maker when it needs to know the
* version number of this file. Since this file is the only unique code in the
* application, the version number of this file is the appropriate one
* to assign to the production attribute "Version". All of the rest
* of the code (the "zeb wrapper") is in the library libBW.a and shouldn't
* change with each measurement. Basically, just leave this function alone
* and everything will work the way it is supposed to.
*****************************************************************************/
{
return (rcsstate);
} /* custom_GetVersion (... */
char*
custom_GetDate ()
/****************************************************************************
* This function is called by dc_maker when it needs to know the
* version date of this file. Since this file is the only unique code in the
* application, the version date of this file is the appropriate one
* to assign to the production attribute "Version_date". All of the rest
* of the code (the "zeb wrapper") is in the library libBW.a and shouldn't
* change with each measurement. Basically, just leave this function alone
* and everything will work the way it is supposed to.
*****************************************************************************/
{
return (rcsdate);
} /* custom_GetDate (... */
void custom_define_hooks ()
{
void parse_commandline_options (void);
/* void lmatherr (struct exception *, int*); */
void create_quicklooks ();
top_hook = parse_commandline_options;
/* matherr_hook = lmatherr; */
/* before_retrieval_hook = change_platforms; */
//before_retrieval_hook = parse_commandline_options;
end_hook = create_quicklooks;
} /* custom_define_hooks (... */
#ifdef BW_STUFF
/***************************************************/
/* Defining a matherr function, so that we can get */
/* sqrt domain errors into the eventLogger */
/***************************************************/
void lmatherr (struct exception *x, int *flag)
{
*flag = 1;
switch (x->type)
{
case DOMAIN:
/* change sqrt to return sqrt(-arg1), not NaN */
if (!strcmp(x->name, "sqrt")) {
x->retval = sqrt (-x->arg1);
msg_ELog (EF_PROBLEM, "Note: SQRT of negative number");
*flag = -1;
} /* if (!strcmp... */
break;
default:
break;
} /* switch (x->type... */
} /* lmatherr (... */
#endif
#define DO_PARSE
#ifdef DO_PARSE
/************************************************************************
* This routine parses the extra flags on the command line, and sets
* some global variables. It then looks to see if the "Plot Only"
* option (-P) was set, and if so, it immediately calls the routine to
* plot the quicklooks (it will not process any data).
************************************************************************/
void parse_commandline_options (void)
{
char **argv;
int argc;
extern int optind, opterr, optopt;
/* extern char *optarg; */
char c;
char datebuf[10]="\0";
char proc_name[30];
/* function defs */
void create_quicklooks (void);
/* Setup for getopt: ignore other commands */
opterr = 0;
platform[0] = '\0';
gLSFITSD = -9999;
/* Get the arguments */
argv = main_GetArgv (&argc);
while ((c = getopt (argc, argv, "n:b:d:l:")) != EOF)
{
switch (c)
{
case 'b':
case 'd':
// Have to pull out first eight characters, which are the date
// This is kludgy, and means we pretty much always have to run
// midnight to midnight, but whatever
strncpy(datebuf,optarg,8);
rundate = atoi (datebuf);
if (rundate < 19000000)
{
rundate += 19000000;
}
break;
case 'n':
#ifdef notdef
sprintf(proc_name, "%s", optarg);
if (strcmp(proc_name, "langley_mfrsr") == 0) {
Platform_List[0] = "mfrsr.b1";
}
#endif
break;
case 'l':
/* LSFITSD (least squars fit standard deviation limit)
*
* Optional command line argument
*
* For some sites, we need to use non-default values
* so make this a command line argument.
*
* This value defaults to 0.006
*
*/
gLSFITSD = atof(optarg);
break;
default: break;
} /* end switch (c) */
} /* end while (c...) */
optind = 1;
if (gLSFITSD <= 0.0) {
gLSFITSD = 0.018; /* Default value */
msg_ELog (EF_INFO, "Default Acceptance Criteria used 0.018\n");
}
} /* parse_commandline_options (... */
#endif
/************************************************************************
* This routine first parses the BW command line to see if the user
* specifically turned the quicklooks off (this is possible with the -N
* command line option). If quicklooks are desired, we then spawn off
* a routine to create the images. This function will need access to
* global variables defined earlier via parse_commandline_options().
************************************************************************/
void create_quicklooks (void)
{
pid_t pid;
int status;
char command[1024], argument[1024];
char setdate[1024],
setplatform[1024], setsite[1024], setfacility[1024];
char **argv, *vapname;
int argc;
void wait();
/*
* Get the arguments
* Note that the only reason I need to do this is to get the name
* of the executable. Typically, I would use main_GetExecName()
* for this task, but if I am here via the top_hook() (i.e., the
* -P option was set), then the variable returned by the
* main_GetExecName function has not been initialized yet...
*/
argv = main_GetArgv (&argc);
vapname = strrchr(argv[0], '/');
if (vapname == NULL) vapname = argv[0];
else vapname = vapname+1;
/* If flag is still true, create the quicklooks */
if (quicklook_flag)
{
msg_ELog (EF_INFO, "Attempting to create quicklooks");
/* Create the child process */
pid = fork ();
if (pid == 0)
{
/*
* If I am here, then I am in the child process.
* Set some environment variables to guide IDL...
*/
/*
* It is important that each environment variable has its own
* memory allocated to it (ie., its own variable). This has
* to do with how putenv() actually works.
*/
sprintf (setdate, "DATE=%d", rundate);
putenv (setdate);
// printf ("DEBUG: DATE = %s\n", setdate);
/* For the ADI version of this VAP, we need to pull out either
* "nimfr" or "mfrsr" from the platform
*
*/
if (strstr(platform, "mfrsr") != (char *) NULL)
{
sprintf (setplatform, "PLATFORM=mfrsr");
}
else
{
sprintf (setplatform, "PLATFORM=nimfr");
}
putenv (setplatform);
// printf ("DEBUG: PLATFORM = %s\n", setplatform);
sprintf (setsite, "SITE=%s", &gSite[0]);
putenv (setsite);
// printf ("DEBUG: SITE = %s\n", setsite);
sprintf (setfacility, "FACILITY=%s", &gFacility[0]);
putenv (setfacility);
// printf ("DEBUG: FACILITY = %s\n", setfacility);
/*
* Create the name of the batch file, as well as the command
* Note that the output from IDL will go to stdout.
*/
sprintf (command, "%s/bin/idl", getenv ("IDL_DIR"));
sprintf (argument, "-em=%s/%s_batch.sav", getenv ("IDL_SAV_DIR"),
vapname);
/*
* Execute the IDL command. If successful, the exit of the
* exec command will terminate the life of the child process
*/
execl (command, "idl", argument, (char *)0);
/*
* This line should never be reached if exec finishes correctly
* We will later trap for this flag to see if everything went
* semi smoothly. Note that we are unable to trap errors that
* might have occurred in IDL when it was running...
*/
exit (127);
} /* end if (pid... */
else
{
/*
* This is the parent process. Wait until the child
* process quits before continuing execution.
*/
wait(&status);
} /* end else */
/* Check the status code from the child process for an error */
if (status == 127)
{
//bw_abort(-1,
msg_ELog(EF_PROBLEM,
"Encountered a problem quicklook generation script/call");
}
} /* end if (quicklook_flag... */
/* If this flag was set, then stop the execution. */
if (plot_only_flag)
{
msg_ELog (EF_INFO, "Program %s completed normally.\n", vapname);
msg_ELog (EF_INFO,
"\n**************************************************************");
exit(0);
} /* end if (plot_only_flag... */
} /* end create_quicklooks(... */
int getStartEndTime (DATA *D,
long *startime,
long *etime)
{
#define SECOFFSET 43200
char srundate[9];
char yyyy[5], mm[3], dd[3];
long timeOffset;
ZebTime midday;
long adjust_to_azimuth(DATA *, long);
/* timeOffset is in hours */
/* This code assumes that lon is east */
if (abs (lon) > 180.0)
{
bw_return (-1, "|longitude| (%f) is greater than 180", lon);
}
timeOffset = (lon == 0.0) ? 0 : (long)(lon / 15.0);
if (timeOffset > 0.0)
{
timeOffset++;
}
sprintf (srundate, "%d", rundate);
strncpy (yyyy, srundate, 4);
strncpy (mm, srundate + 4, 2);
strncpy (dd, srundate + 6, 2);
yyyy[4] = '\0';
mm[2] = '\0';
dd[2] = '\0';
TC_ZtAssemble (&midday, atoi (yyyy), atoi (mm), atoi (dd), 12, 0, 0, 0);
// Uncomment out the following four lines to remove azimuth correction
//midday.zt_Sec -= (timeOffset *3600);
//*startime=midday.zt_Sec-SECOFFSET;
//*etime=midday.zt_Sec+SECOFFSET;
//return(TRUE);
/* convert midday GMT to midday local time */
long noon = midday.zt_Sec - (timeOffset *3600);
long ta,tb;
double aza,azb;
double get_azimuth(long);
// find midnight before and midnight after midday, and that's our start
// to end time. Jump 10 hours, and then go minute by minute. The test
// for crossing midnight is az(ta) > 180 && az(tb) < 180, for tb>ta
ta=noon-10*3600;
azb=get_azimuth(ta+60);
aza=get_azimuth(ta);
// prev midnight
while (! (aza > 180. && azb <= 180) ) {
ta -= 60;
azb=aza;
aza=get_azimuth(ta);
}
*startime = ta+30; // tb is in , ta is out
// following midnight
tb=noon+10*3600;
azb=get_azimuth(tb);
aza=get_azimuth(tb-60);
while (! (aza > 180. && azb <= 180) ) {
tb += 60;
aza=azb;
azb=get_azimuth(tb);
}
*etime = tb-30; // ta is in , tb is out
return (TRUE);
} /* end getStartEndTime */
// Function to calculate azimuth at a given time (for lat and lon)
double get_azimuth(long et) {
ZebTime zt;
UItime dt;
double efactor, el, ap_ra, ap_dec, refrac, distance;
double az;
double dayfrac;
extern int solarposition();
int yymmdd_to_juldate();
int hh, min, ss, jday, year;
zt.zt_Sec = et;
zt.zt_MicroSec = 0;
// Calculate azimuths, which is a whole big proces
TC_ZtToUI(&zt, &dt);
jday = yymmdd_to_juldate(dt.ds_yymmdd);
year = dt.ds_yymmdd / 10000;
if (year < 1900) year += 1900; /* y2k or not */
hh = dt.ds_hhmmss / 10000;
min = (dt.ds_hhmmss /100 ) % 100;
ss = dt.ds_hhmmss % 100;
dayfrac = (double) (hh*3600 + min*60 + ss) / (24.0*3600.0);
dayfrac += (double) jday;
double arg;
if (year % 4 == 0)
arg = (double) (jday - 1)/366.0*2.0*M_PI;
else
arg = (double) (jday - 1)/365.0*2.0*M_PI;
efactor = 1.00011 + .034221*cos(arg) + .00128*sin(arg) +
.000719*cos(2*arg) + 7.7e-5*sin(2*arg);
solarposition(year, 0, dayfrac, (double) 0.0, lat, lon,
&ap_ra, &ap_dec, &el, &refrac, &az, &distance);
return(az);
}
void checkQCFlag (DATA *D)
{
int o, f, t;
int qc_value;
/*
* loop through all obs and samples and if the qc flag is set
* set the corresponding data value to MISSING.
*/
// printf ("DEBUG: nSamples = %d\n", D->nSamples[B1][0][0][0]);
// printf ("DEBUG: NCHANNELS = %d\n", NCHANNELS);
for (o = 0; o < D->nObs[B1]; o++)
{
for (t = 0; t < D->nSamples[B1][o][0][0]; t++)
{
/* loop though the chn fields */
for (f = 0; f < (NCHANNELS - 1); f++)
{
switch(f)
{
case 0:
/* p o [0] f t n */
qc_value = (int) D->BWdata[B1][o][0][QC_CHN1][t][0];
memcpy (&qc_value,
(int *) &D->BWdata[B1][o][0][QC_CHN1][t][0],
sizeof(int));
if (qc_value != 0)
{
D->BWdata[B1][o][0][IN_CHN1][t][0] = MISSING;
}
break;
case 1:
memcpy (&qc_value,
(int *) &D->BWdata[B1][o][0][QC_CHN2][t][0],
sizeof(int));
if (qc_value != 0)
{
D->BWdata[B1][o][0][IN_CHN2][t][0] = MISSING;
}
case 2:
memcpy (&qc_value,
(int *) &D->BWdata[B1][o][0][QC_CHN3][t][0],
sizeof(int));
if (qc_value != 0)
{
D->BWdata[B1][o][0][IN_CHN3][t][0] = MISSING;
}
break;
case 4:
memcpy (&qc_value,
(int *) &D->BWdata[B1][o][0][QC_CHN4][t][0],
sizeof(int));
if (qc_value != 0)
{
D->BWdata[B1][o][0][IN_CHN4][t][0] = MISSING;
}
break;
case 5:
memcpy (&qc_value,
(int *) &D->BWdata[B1][o][0][QC_CHN5][t][0],
sizeof(int));
if (qc_value != 0)
{
D->BWdata[B1][o][0][IN_CHN5][t][0] = MISSING;
}
break;
case 6:
memcpy (&qc_value,
(int *) &D->BWdata[B1][o][0][QC_CHN6][t][0],
sizeof(int));
if (qc_value != 0)
{
D->BWdata[B1][o][0][IN_CHN6][t][0] = MISSING;
}
break;
}
} /* end for (f... */
/* check the broadband field */
memcpy (&qc_value,
(int *) &D->BWdata[B1][o][0][QC_BROADBAND][t][0],
sizeof(int));
if (qc_value != 0)
{
D->BWdata[B1][o][0][IN_BROADBAND][t][0] = MISSING;
}
} /* end for (t... */
} /* end for (o... */
} /* end checkQCFlag (... */
void getWavelengthAttr (DATA *D)
{
char *foo;
int r;
/*
* 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++)
{
foo = (char *)bw_GetAttrValue (D,
B1,
0,
IN_CHN1 + r,
"actual_wavelength",
NULL,
NULL);
if (foo != NULL)
{
channels[r + 1] = atof (foo);
// printf ("DEBUG: r = %d ... actual_wavelength = %f\n", r, atof(foo));
}
} /* end for (r... */
channels[0] = 0; /* This is the broadband channel */
/* return (TRUE); */
} /* end getWavelengthAttr (... */
int writeWavelengthAttr (DATA *newD)
{
char chns[16];
int o, r;
for (o = 0; o < newD->nObs[LANGLEY]; o++)
{
for (r = 0; 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]);
}
/* barnard */
bw_SetAttrValue (newD, LANGLEY, o, OD + r, DCT_String,
1, "actual_wavelength", (void *)chns);
bw_SetAttrValue (newD, LANGLEY, o, SC + r, DCT_String,
1, "actual_wavelength", (void *)chns);
bw_SetAttrValue (newD, LANGLEY, o, ERR_FIT + r, DCT_String,
1, "actual_wavelength", (void *)chns);
bw_SetAttrValue (newD, LANGLEY, o, ERR_SLOPE + r, DCT_String,
1, "actual_wavelength", (void *)chns);
bw_SetAttrValue (newD, LANGLEY, o, NPTS + r, DCT_String,
1, "actual_wavelength", (void *)chns);
bw_SetAttrValue (newD, LANGLEY, o, GOOD_FRACT + r, DCT_String,
1, "actual_wavelength", (void *)chns);
bw_SetAttrValue (newD, LANGLEY, o, BAD_FLAG + r, DCT_String,
1, "actual_wavelength", (void *)chns);
/* michalsky */
bw_SetAttrValue (newD, LANGLEY, o, OD + NCHANNELS + r, DCT_String,
1, "actual_wavelength", (void *)chns);
bw_SetAttrValue (newD, LANGLEY, o, SC + NCHANNELS + r, DCT_String,
1, "actual_wavelength", (void *)chns);
bw_SetAttrValue (newD, LANGLEY, o, SC_NORM + r, DCT_String,
1, "actual_wavelength", (void *)chns);
bw_SetAttrValue (newD, LANGLEY, o, SD + r, DCT_String,
1, "actual_wavelength", (void *)chns);
bw_SetAttrValue (newD, LANGLEY, o, NPTS + NCHANNELS + r, DCT_String,
1, "actual_wavelength", (void *)chns);
bw_SetAttrValue (newD, LANGLEY, o, GOOD_FRACT + NCHANNELS + r,
DCT_String, 1, "actual_wavelength", (void *)chns);
bw_SetAttrValue (newD, LANGLEY, o, BAD_FLAG + NCHANNELS + r, DCT_String,
1, "actual_wavelength", (void *)chns);
} /* end for (r... */
} /* end for (o... */
for (o = 0; o < newD->nObs[PLOT]; o++)
{
for (r = 0; 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]);
}
/* barnard */
bw_SetAttrValue (newD, PLOT, o, LN_I + r, DCT_String,
1, "actual_wavelength", (void *)chns);
bw_SetAttrValue (newD, PLOT, o, REJECTED + r, DCT_String,
1, "actual_wavelength", (void *)chns);
/* michalsky */
bw_SetAttrValue (newD, PLOT, o, LN_I + NCHANNELS + r, DCT_String,
1, "actual_wavelength", (void *)chns);
bw_SetAttrValue (newD, PLOT, o, REJECTED + NCHANNELS + r, DCT_String,
1, "actual_wavelength", (void *)chns);
} /* end for (r... */
} /* end for (o... */
return (TRUE);
} /* writeWavelengthAttr (... */
DATA*
readFromExternal (char* externalInFile)
/*****************************************************************************
* Input parameters:
*
* Returned value:
*
* Called from:
*
* Functions called:
*
* Description:
* This function is called only for data ingest. In that case, the
* programmer must write the code in this function that will read in
* the desired data (from the argument externalFile) and then pack the
* data into the DATA structure d, which is returned. Code should not
* entered here for measurement and data delivery applications.
*
*
*****************************************************************************/
{
/* Local variables */
DATA *d = NULL;
return (d);
}
#undef LANGLEY_C
/* langley.c */
src/langley.h 0 → 100644
View file @ 02287497
/*******************************************************************************
* COPYRIGHT (C) 2001 Battelle Memorial Institute.
* All Rights Reserved. (Now, with this out of the way...)
*
* RCS INFORMATION:
* $RCSfile: langley.h,v $
* $Revision: 1.15 $
* $Author: koontz $
* $Locker: $
* $Date: 2010/01/22 20:45:31 $
* $State: Exp $
* $Name: $
* $Id: langley.h,v 1.15 2010/01/22 20:45:31 koontz Exp $
*
* LIBRARY DEPENDENCIES:
*
* DESCRIPTION:
* <Discuss overall design/purpose of class.>
*******************************************************************************/
#ifndef LANGLEY_H
#define LANGLEY_H
/* We'll be doing a lot of dynamic allocation, so let's make it easier */
#define CALLOC(n,t) (t*)calloc(n, sizeof(t))
#define REALLOC(p,n,t) (t*) realloc((char *) p, (n)*sizeof(t));
#define MISSING -9999.0
/*
* Really should have number of narrow and number of broad band channels,
* and have a narrow_idx[*] and a broad_idx[*] and work off of the idxs.
* But for now we will just assume that the fist chn is the broad and
* narrow starts at 1
*/
#define NCHANNELS 7 /* number of mfrsr channels 1 broadband, 6 narrow band */
#define KEY_CHN 1 /* this is the 500nm channel */
#define NSLABS 4 /* max number of langleys per day, should only get 2 */
#define PLOTDATA /* if set generate the plot netcdf file. */
#define LOW_AM 2.0
#define HIGH_AM 6.05
#define LSFITSD 0.006
#define OUTLIMIT 1.5
#define FRACPTS 0.3
#define CLOUDSLOP 0.015
#define DEFAULT_PLATFORM "mfrsr"
#define DEFAULT_SITE "sgp"
#define DEFAULT_FACILITY "C1"
#define UNCAL_DATALEVEL "a0"
#define UNCAL_NAME "uncalibrated"
enum platform_type
{
Input_Platform,
Output_Platform
}; /* enum platform_type */
enum input_platforms
{
B1_M = 0, /* platform #0 in retriever file is the mfrsr b1 datastream */
B1_N = 1, /* platform #0 in retriever file is the nimfr b1 datastream */
}; /* enum input_platforms */
int B1;
#define NPRESPLATS 4
#define PRES_CONF_FILE "langley_pres.conf"
enum input_fields
{
IN_BROADBAND = 0,
IN_CHN1,
IN_CHN2,
IN_CHN3,
IN_CHN4,
IN_CHN5,
IN_CHN6,
QC_BROADBAND,
QC_CHN1,
QC_CHN2,
QC_CHN3,
QC_CHN4,
QC_CHN5,
QC_CHN6
}; /* enum input_fields */
enum input_pres_flds
{
PRES
}; /* enum input_pres_flds */
enum output_fields
{
OD = 0, /* barnard then michalsky langley */
SC = 14, /* barnard then michalsky langley */
SC_NORM = 28, /* michalsky langley */
ERR_FIT = 35, /* barnard langley */
SD = 42, /* michalsky langley */
ERR_SLOPE = 49, /* barnard langley */
NPTS = 56, /* barnard then michalsky langley */
GOOD_FRACT = 70, /* barnard then michalsky langley */
BAD_FLAG = 84, /* barnard then michalsky langley */
AM = 0, /* barnard then michalsky plot */
LN_I = 2, /* barnard then michalsky plot */
BSDIST = 98, /* barnard sun earth distance */
MSDIST = 99, /* michalsky sun earth distance */
REJECTED = 16 /* barnard then michalsky plot */
}; /* enum input_fields */
#define NOUTFLDS 100
#define NOUTPLTFLDS 30
#define NOUTPLATS 2
enum output_platforms
{
LANGLEY,
PLOT
}; /* enum output_platforms */
#ifdef LANGLEY_C
char altPlatform[9] = "Altitude";
float channels[NCHANNELS] = {413., 499., 608., 665., 859., 939., 0};
double lat, lon, alt;
/* double presfraction; */
#else
extern char altPlatform[9];
extern int channels[NCHANNELS];
extern double lat, lon, alt;
/* extern double presfraction; */
#endif /* LANGLEY_C */
float gLSFITSD;
/****** function prototypes ******/
#endif /* LANGLEY_H */
/* LANGLEY_H */
src/langley.pro 0 → 100644
View file @ 02287497
;-----------------------------------------------------------------------------
; COPYRIGHT (C) 2001 Battelle Memorial Institute.
; All Rights Reserved. (Now, with this out of the way...)
;
; RCS INFORMATION:
; $RCSfile: langley.pro,v $
; $Revision: 1.17 $
; $Author: koontz $
; $Locker: $
; $Date: 2012-01-09 16:02:40 $
; $State: Exp $
; $Name: not supported by cvs2svn $
; $Id: langley.pro,v 1.17 2012-01-09 16:02:40 koontz Exp $
;
; FUNCTIONS IN THIS FILE:
;
; ABSTRACT:
; This script does the Langley plots and marks the rejected points
; using data from the sgplangleyXX.c1 and sgplangplotXX.c1 platforms.
; Possibly the worst written IDL script in the history of the universe.
;
; The main procedure is "dolangley", which generates a Langley plot from
; sgplangleyC1.c1 ("analysis") and sgplangplotC1.c1 ("plot") data files.
; It loops over all the plot files in the "plotdir" directory; if no
; corresponding analysis file exists in "analdir", the code will exit.
;
; Calling Sequence:
; langley, [<channel>]
;
; Arguments
; channel - the MFRSR channel which to plot. Values given by integers
; 0-6; "0" indicates the broadband channel, "1" is 415 nm,
; etc. If channel < 0 or is missing, all channels will be
; plotted.
;
; Keywords
; date - a string that contains a date or part of one, used to
; restrict which plot files we loop over. As the name
; suggests, this is usually used to plot just one day or month
; of data. If not set, all the files in "plotdir" will be
; looped over.
;
; to - If set, do plots for all the channels between the <channel>
; argument and the value of this keyword.
;
; site - The site for the data (SGP, NSA, TWP)
;
; facility - The facility for the data (B1,B4, B5, B6, C1, etc.)
;
; algorithm - sets the algorithm to use. Can be set to barnard or
; michalksy. Default is barnard.
;
; uncalibrated - set if the plot will use the uncalibrated data.
;
; onefile - If set, all the postscript plots (i.e. for all channels)
; will be sent to one file for each time. This lets you
; use psnup or someting similar to put multiple channel
; plots for the same day on one sheet of paper. Default is
; to write each Langley plot to a seperate postscript file.
; Has no effect unless !D.name is set to 'PS' (i.e. via
; set_plot, 'PS' ).
;
; noplot - If set, no plots are created. Probably only makes sense
; if you also set the /stats keyword, but, hey, go nuts.
;
; stats - If set, generates statistics on number of accepted plots
; for each channel as well number of rejected plots for each
; rejection criterion. Prints a table to stdout once all
; plots are finished; with /noplot, only the statistics will
; be done.
;
; quiet - If set, don't print status updates or names of output
; files; will still print final stats if /stats is set.
;
; 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
; x.
;
; no_catch - for debuging. set to turn off error catching
;
; zlog - set if messages are to be sent to the zebra logger.
;
; Examples
; Do Langley plots for one (prompted) channel for all the files
; currently residing in "plotdir":
;
; IDL> langley
;
; Do Langley plots for channels 1, 2, 3, and 4, for files with "9612"
; somewhere in their name, with all plots going to one big file in
; the end and a table of statistics sent to stdout after all the
; plotting is done:
;
; IDL> set_plot, 'PS'
; IDL> dolangley, 1, to=4, date="9612", /onefile, /stats
;-
pro langley, ch, to = to, date = td, $
platform = platform, site = site, facility = facility, $
algorithm = algorithm, uncalibrated = uncalibrated, $
onefile = onefile, noplot = noplot, stats = stats, $
quiet = quiet, plot_to = plot_to, no_catch = no_catch, $
zlog = zlog
;
; force this rascal to use the Z buffer (so quicklooks will
; work in production)
;
set_plot, 'Z'
; if zlog is set then open a file to print msgs, else print to stdio
if (keyword_set (zlog)) then begin
zlog = 2
msgfile = getenv ('TMP_DATA') + '/' + site + platform + $
'Langley' + facility + '.log'
; set up command to print to the zebra logger
; spawn uses the startup env so we need to set the ZEBRA_SOCKET var
; to what we are currently using.
zlog_command = 'setenv ZEBRA_SOCKET ' + getenv ('ZEBRA_SOCKET') + $
'; zlog -p langley_ql -P -f ' + msgfile
openw, zlog, msgfile, error = ferr
if (ferr ne 0) then begin
printf, ­2, !ERR_STRING
zlog = -1
endif
endif else zlog = -1
; Capture some info about the arguments we got
printf, zlog, ' Langley plot starting execution'
printf, zlog, ' -------------------------------'
printf, zlog, ' Site = ', site
printf, zlog, ' Platform = ', platform
printf, zlog, ' Facility = ', facility
printf, zlog, ' Algorithm = ', algorithm
printf, zlog, ' Date = ', td
printf, zlog, ' Plot_to = ', plot_to
; debug info for DMF operator
if (zlog ne -1) then begin
print, ' Langley plot starting execution'
print, ' -------------------------------'
print, ' Site = ', site
print, ' Platform = ', platform
print, ' Facility = ', facility
print, ' Algorithm = ', algorithm
print, ' Date = ', td
print, ' Plot_to = ', plot_to
close, zlog, /force
openw, zlog, msgfile, error = ferr, /append
printf, zlog, ' reopened log file at A'
print, ' reopened log file at A'
endif
Dname = !D.name
plot_to = "png"
printf, zlog, ' Forcing Langley to use PNG plot type'
if (zlog ne -1) then begin
print, ' Forcing Langley to use PNG plot type'
close, zlog, /force
openw, zlog, msgfile, error = ferr, /append
endif
; Catch error and exit
if (not keyword_set(no_catch)) then begin
catch, error
if (error ne 0) then begin
printf, zlog, ' Langley plot has Aborted'
printf, zlog, ' '+!ERR_STRING
printf, zlog, ' Langley plot got no_catch error ???'
if (zlog ne -1) then begin
close, zlog
print, ' .... zlog_command is ', zlog_command
spawn, zlog_command
endif
return
endif
endif
; Revision number of this file. It shows up in tiny print on the Langley
; plots themselves.
;idlVersion = '8.0'
idlVersion = '8.0'
qlDate = '$Date: 2012-01-09 16:02:40 $'
first = strpos (qlDate, ' ') + 1
last = strpos (qlDate, ' ', /reverse_search)
qlDate = strmid (qlDate, first, last - first)
qlVersion = '$Revision: 1.17 $'
first = strpos (qlVersion, ' ') + 1
last = strpos (qlVersion, ' ', /reverse_search)
qlVersion = strmid (qlVersion, first, last - first)
vapVersion = '$State: Exp $'
first = strpos (vapVersion, ' ') + 1
last = strpos (vapVersion, ' ', /reverse_search)
vapVersion = strmid (vapVersion, first, last - first)
printf, zlog, ' qlDate = ', qlDate
printf, zlog, ' qlVersion = ', qlVersion
printf, zlog, ' vapVersion = ', vapVersion
if (zlog ne -1) then begin
close, zlog, /force
openw, zlog, msgfile, error = ferr, /append
endif
; Input netCDF directories
; These need to be set correctly by you to where you keep the plot and
; analysis data files. Note that you need the "/" at the end.
datastream_home = getenv ('DATASTREAM_DATA')
if (not keyword_set (platform)) then platform = "mfrsr"
if (not keyword_set (site)) then site = "sgp"
if (not keyword_set (facility)) then facility = "C1"
if (not keyword_set (algorithm)) then algorithm = "barnard"
if (keyword_set (uncalibrated)) then begin
uncal = "uncalibrated"
endif else begin
uncal = ""
endelse
anal_ds = site + uncal + platform + 'langley' + facility + '.c1'
anal_ql = site + uncal + platform + 'langley' + algorithm + facility + '.c1'
plot_ds = site + uncal + platform + 'langplot' + facility + '.c1'
plot_dir = datastream_home + '/' + site + '/' + plot_ds + '/'
anal_dir = datastream_home + '/' + site + '/' + anal_ds + '/'
; Initialize some tags
Newfile = 1
MICHALSKY = 1
BARNARD = 2
algorithm = strlowcase (algorithm)
if (algorithm eq "michalsky") then begin
algor = MICHALSKY
endif else begin
if (algorithm eq "barnard") then begin
algor = BARNARD
endif else begin
printf, zlog, "Unknown algorithm - " + algorithm
return
endelse
endelse
;; set the colors
;; set up color table
@color_syms.include
; setup the output type
Dname = !D.name
; if (not keyword_set (plot_to)) then begin
plot_to = "png"
printf, zlog, ' Forcing Langley to use PNG plot type'
; endif
; if (strlowcase (plot_to) eq "png") then begin
printf, zlog, ' Langley plot_to keyword is PNG '
print, Dname
set_plot, 'Z'
print, Dname
device, set_resolution = [1024, 786]
f_tau = "!7s!X!N"
f_sigma = "!7r!X!N"
Font = -1
Bcolor = d_white
Scolor = d_black
; endif
; endif else begin
; if (strlowcase (plot_to) eq "ps") then begin
; printf, zlog, ' Langley plot_to keyword is PS '
; set_plot, 'PS'
; ; ps device is set when file name is set
; f_tau = "!Mt!X!N"
; f_sigma = "!Ms!X!N"
; Font = 0
; Bcolor = d_white
; Scolor = d_black
; endif else begin
; if (strlowcase (plot_to) eq "x") then begin
; set_plot, 'Z'
; set_plot, 'X'
; xsize = 1024
; ysize = 768
; f_tau = "!7s!X!N"
; f_sigma = "!7r!X!N"
; Font = -1
; Bcolor = d_white
; Scolor=d_black
; endif
; endelse
; endelse
; These shouldn't be hardcoded; we should be able to pull it out from
; the data itself. Oh well.
NCHNS = 7
chname = ['broadband', '415 nm', '499 nm', '608 nm', '664 nm', $
'860 nm', '938 nm']
if (algor eq MICHALSKY) then begin
rejwhy = ['no points to regress', 'all points filtered', $
'???']
endif else begin
rejwhy = ['error in fit > .01', 'solar constant out of range', $
'good fraction < .75', 'chn good, marked bad by key chn']
rejwhy_415 = ['error in fit > .015', 'solar constant out of range', $
'good fraction < .75', 'chn good, marked bad by key chn']
endelse
; Make our stats array if /stats is set
if (n_elements(stats) gt 0) then begin
nplots = make_array(NCHNS, /integer, value=0)
badflag_counter = make_array(NCHNS, 4, /integer, value=0)
endif
; Funky initializations if /all is set; we have to set ch to something to
; foil the prompting that happens right after this.
if (n_elements(all) gt 0) then begin
ch = -1
td = ''
endif
; Prompt a channel to display, if not given
if (n_elements(ch) eq 0 ) then ch = -1
; Negative channel means do them all
if (ch le 0) then begin
ch = 0
endch = 6
endif else begin
endch = ch
endelse
; Set the end channel if /to is set
if (n_elements (to) ne 0) then endch = to
; Set the selection string to nothing, if not given
if (n_elements (td) eq 0) then td = ''
td = strtrim (td, 2)
; This scans through all available files, and uses the selection string
lookfor = plot_dir + '*' + td + '*.cdf'
print, "looking for files: " + lookfor
ls = findfile (plot_dir + '*' + td + '*.cdf', count = nfoo)
if (n_elements (quiet) eq 0) then printf, zlog, $
'Running over ' + string(nfoo) + ' plots...'
for i=0, nfoo-1 do begin
;; Get this plot file
fnew=ls(i)
;; We now have the langplot file name - so open it up
fin = ncdf_open(fnew)
;; read in the times
BaseTime = 0L
ncdf_varget, fin, 'base_time', BaseTime
ncdf_varget, fin, 'time_offset', TimeOffset
;; convert to julian time - we'll need this later when we try to
;; track down the right analysis data for this plot.
t_zt2julian, BaseTime, TimeOffSet, jday
nslab = n_elements(jday)
;; get our YYMMDD.hhmmss time stamp, for headers and file names
t_zt2hhmmss, BaseTime, TimeOffset, dateYYMMDD, timeHHMMSS
Tstamp = strtrim (dateYYMMDD(0), 2) + '.' + strtrim (timeHHMMSS(0), 2)
if (n_elements (quiet) eq 0) then printf, zlog, $
'==> ' + Tstamp + ' <=='
;; Now read in the important data
lograd = dblarr (n_elements (TimeOffset), NCHNS)
rejected = lograd
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 (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
lograd[*, j] = foonew
foonew = 0
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
rejected[*, j] = foonew
foonew = 0
endfor
lograd = transpose (lograd)
rejected = transpose (rejected)
; get input files version info
ncdf_attget, fin, 'process_version', finVersion, /global
finVersion = string (finVersion)
first = strpos (finVersion, ' ') + 1
last = strpos (finVersion, ' ', /reverse_search)
finVersion = strmid (finVersion, first, last - first)
ncdf_attget, fin, 'history', finDate, /global
finDate = string (finDate)
first = strpos (finDate, ' at ') + 4
last = strpos (finDate, ', using ', /reverse_search)
finDate = strmid (finDate, first, last - first)
;; close the netcdf file
breakout:
ncdf_close, fin
;; Now grab the analysis file.
dateYYMMDD = strtrim (dateYYMMDD, 2)
ymdhms2systime, strmid (dateYYMMDD, 0, 4), $
strmid (dateYYMMDD, 4, 2), $
strmid (dateYYMMDD, 6, 2), 0, 0, 0, mdate
sdate = mdate - (24 * 60 * 60)
nbar = 0
for d = 0, 2 do begin
t_zt2hhmmss, sdate + (d * 24 * 60 * 60), 0, newYYMMDD
lfiles = findfile (anal_dir + '*' + strtrim (newYYMMDD, 2) + '*.cdf', $
count = nbar)
k = 0
;; scan down all the analysis files for this day, and find one where
;; a sample time is in the middle of the period covered by the plot file
;; That's the analysis sample we want.
for j=0, nbar-1 do begin
flang = ncdf_open(lfiles(j))
BT = 0L
ncdf_varget, flang, 'base_time', BT
ncdf_varget, flang, 'time_offset', TO
t_zt2julian, BT, TO, jules
for k=0, n_elements(jules)-1 do begin
if (jules(k) le jday(nslab-1) and $
jules(k) ge jday(0)) then goto, read_flag
endfor
;; This means this analysis file is not the one we want, so close it
ncdf_close, flang
endfor
endfor
;; if we're here, that means there is no analysis file for this day
;; this is probably a serious problem
printf, zlog, 'No analysis file here? Huh?'
exit
;; k should be the appropriate index for the open flang file
read_flag:
t_zt2hhmmss, BT, TO, dyymmdd, thhmmss
tau = dblarr (n_elements (TO), NCHNS)
I0 = tau
sd = tau
sigs = tau
npoints = tau
frac = tau
bad = tau
if (algor eq MICHALSKY) then begin
scvar = 'michalsky_solar_constant_sdist'
sdvar = 'michalsky_standard_deviation'
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
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
for j=1, NCHNS - 1 do begin
ncdf_varget, flang, algorithm + '_optical_depth_filter' + $
strtrim (j, 2), foo
tau[*, j] = foo
foo = 0
ncdf_varget, flang, scvar + '_filter' + strtrim (j, 2), foo
I0[*, j] = foo
foo = 0
ncdf_varget, flang, sdvar + '_filter' + strtrim (j, 2), foo
sd[*, j] = foo
foo = 0
ncdf_varget, flang, algorithm + '_number_points_filter'+ $
strtrim (j, 2), foo
npoints[*, j] = foo
foo = 0
ncdf_varget, flang, algorithm + '_good_fraction_filter' + $
strtrim (j, 2), foo
frac[*, j] = foo
foo = 0
ncdf_varget, flang, algorithm + '_badflag_filter' + $
strtrim (j, 2), foo
bad[*, j] = foo
foo = 0
if (algor eq BARNARD) then begin
ncdf_varget, flang, 'barnard_error_slope_filter' + strtrim (j, 2), foo
sigs[*, j] = foo
foo = 0
endif
endfor
tau = transpose (tau)
I0 = transpose (I0)
sd = transpose (sd)
sigs = transpose (sigs)
npoints = transpose (npoints)
frac = transpose (frac)
bad = transpose (bad)
; get input files version info
ncdf_attget, fin, 'process_version', flangVersion, /global
flangVersion = string (flangVersion)
first = strpos (flangVersion, ' ') + 1
last = strpos (flangVersion, ' ', /reverse_search)
flangVersion = strmid (flangVersion, first, last - first)
ncdf_attget, fin, 'history', flangDate, /global
flangDate = string (flangDate)
first = strpos (flangDate, ' at ') + 4
last = strpos (flangDate, ', using ', /reverse_search)
flangDate = strmid (flangDate, first, (last - first))
;; close 'er up
ncdf_close, flang
;; Now loop over all the channels we have selected
for channel=ch, endch do begin
if (n_elements(quiet) eq 0) then printf, zlog, $
'Plotting ' + chname(channel) + ' channel'
;; Restrict to the channel we are currently looping over.
A = airmass
lnI = lograd(channel, *)
rej = rejected(channel, *)
;; convert analysis products into something we understand
slope = - tau(channel, k)
intercept = alog(I0(channel, k))
sdp = sd(channel, k)
if (algor eq MICHALSKY) then begin
error_slope = -9999.0
endif else begin
error_slope = sigs[channel, k]
endelse
num = npoints(channel, k)
good_fraction = frac(channel, k)
badflag = bad(channel, k)
;; Here is our statistical bookkeeping
if (n_elements(stats) gt 0) then begin
nplots(channel) = nplots(channel) + 1
badflag_counter(channel, fix(badflag)) = $
badflag_counter(channel, fix(badflag)) + 1
endif
;; skip all the actual plotting stuff if you've set /noplot
if (n_elements(noplot) gt 0) then goto, skip_plot
;; Set the window number for the plot
winNum = 0
;; Title of the plot
Ptitle = 'Langley plot for ' + chname(channel) + '!C' + Tstamp
if (!D.name eq 'PS') then begin
get_quicklook_dir, dateYYMMDD, uncal + platform + 'langley', site, $
facility, 'c1', /create, $
quicklook_dir = ql_dir, errornum = errnum
ql_dir = ql_dir + '/'
;; This conditional tells us whether to open a new file or not
;; If /onefile is not set, you always open a new file;; otherwise
;; you only do it once (Newfile = 1 at the start).
if (n_elements(onefile) eq 0 or Newfile gt 0) then begin
;; we also add an "all" to the output filename when we set
;; /onefile so it looks different than individual plot files.
if (n_elements(onefile) ne 0) then begin
fname = ql_dir + anal_ql + ".ch" + $
string(format='(I1.1)', channel) + '-' + $
string(f='(I1.1)',endch) + '.' + Tstamp + ".ps"
endif else begin
fname = ql_dir + anal_ql + ".ch" + $
string(format='(I1.1)', channel) + '.' + $
Tstamp + ".ps"
endelse
if (n_elements(quiet) eq 0) then printf, zlog, $
"Output file is " + fname
;; setup the postscript file
device, filename=fname, /times, /inches, /color, /landscape
endif
;; Reset newfile if onefile is set and we are on the last channel
;; This will cause a new plot file the next time through
if (n_elements(onefile) ne 0 and channel ne endch) then begin
Newfile = 0
endif else begin
Newfile = 1
endelse
endif
;; open a window
if (!d.name eq 'X') then window, winNum, retain=2, $
xsize=xsize, ysize=ysize
!p.multi= 0
!p.region=[0,0.05,1.0,0.95]
!p.font = Font
!p.color = Scolor
!p.background = Bcolor
;; This is clunky - 'foo' gives indeces of accepted points, 'bar' has
;; the indeces of rejected points. I could have named them better.
foo = where(rej eq 0)
bar = where(rej ne 0)
;; finally, do the plots
plot, A, lnI, color=Scolor, $
min_value = -9000, $
xrange = [2, 6.05], $
ytitle = "log(irradiance)", xtitle = "Airmasses", $
title=Ptitle,charsize=1.5, xstyle=1, ystyle=1, /nodata
;; accepted points are green plus signs
if (foo(0) gt -1) then begin
oplot, A(where(rej eq 0)), lnI(where(rej eq 0)), $
color=d_green, psym=1
endif
;; rejected points are red crosses
if (bar(0) gt -1) then begin
oplot, A(where(rej ne 0)), lnI(where(rej ne 0)), $
color=d_red, psym=7
endif
;; plot the regression line, with the two sigma lines in cyan
amin = min(A)
amax = max(A)
oplot, [amin, amax], [intercept+slope*amin, intercept+slope*amax], $
color = Scolor
oplot, [amin, amax], $
[intercept+slope*amin-2*sdp, $
intercept+slope*amax-2*sdp], $
color = d_cyan, linestyle=2
oplot, [amin, amax], [intercept+slope*amin+2*sdp, $
intercept+slope*amax+2*sdp], $
color = d_cyan, linestyle=2
;; Say whether the whole plot is rejected (in red) or accepted
;; (in green)
if (badflag eq 0.0) then begin
xyouts, 0.05, .95, "ACCEPTED", color = d_green, charsize=2.0, /normal
endif else begin
xyouts, 0.05, .95, $
'REJECTED(' + string(format='(I1.1)', badflag) + ')', $
color = d_red, charsize=2.0, /normal
if (algor eq MICHALSKY) then begin
xyouts, 0.05, .93, rejwhy(badflag-1), color = d_red, /normal
endif else begin
;; slightly different rejection(1) criteria for the 415 line,
;; so we have to switch on the channel here.
if (channel eq 1) then begin
xyouts, 0.05, .93, rejwhy_415(badflag-1), color = d_red, /normal
endif else begin
xyouts, 0.05, .93, rejwhy(badflag-1), color = d_red, /normal
endelse
endelse
endelse
;; Now write out the other plot information: optical depth tau,
;; error in tau, error in fit, solar constant I0, # of (good)
;; points used in final regression, percentage of total points
;; that are good
;; get the y limits on the plot, so we can put our text in the
;; right place
idxNotMissing = where (lnI gt -9999)
if (idxNotMissing[0] eq -1) then begin
imin = 0
imax = 1
endif else begin
imin = min(lnI(idxNotMissing))
; imin = min(lnI(where(lnI gt -9999)))
imax = max(lnI)
endelse
!p.charsize = 1.5
xyouts, 0.89, 0.95, f_tau + ":", $
alignment = 1, /normal
xyouts, 0.90, 0.95, string(format='(F6.4)', -slope), $
alignment = 0, /normal
xyouts, 0.89, 0.925, f_sigma + "!D" + f_tau + "!X!N:", $
alignment = 1, /normal
xyouts, 0.90, 0.925, string(format='(F6.4)', error_slope), $
alignment = 0, /normal
!p.charsize = 1.2
if (algor eq MICHALSKY) then begin
xyouts, 5.9, imin+(imax-imin)*.95,"standard deviation = " + $
string(format='(F8.5)', sdp), alignment=1
endif else begin
xyouts, 0.85, 0.85, "error in fit:", alignment = 1, /normal
xyouts, 0.86, 0.85, string(format='(F7.5)', sdp), $
alignment = 0, /normal
endelse
xyouts, 0.85, 0.82, "I!D0!N:", alignment=1, /normal
xyouts, 0.86, 0.82, string(format='(F7.5)', I0(channel, k)), $
alignment=0, /normal
xyouts, 0.85, 0.79,"fraction good:", alignment = 1, /normal
xyouts, 0.86, 0.79, string(format='(F7.5)', good_fraction), $
alignment = 0, /normal
xyouts, 0.85, 0.73, "# of points:", alignment = 1, /normal
xyouts, 0.86, 0.73, string(format='(I4)', num), $
alignment = 0, /normal
; print out algorithm info
!p.charsize = 1.0
xyouts, 0.91, 0.07, "Algorithm:", alignment = 1, /normal
xyouts, 0.92, 0.07, algorithm, alignment = 0, /normal
xyouts, 0.91, 0.05, "Location:", alignment = 1, /normal
xyouts, 0.92, 0.05, site + ' ' + facility, alignment = 0, /normal
xyouts, 0.91, 0.03, "DataType:", alignment = 1, /normal
if (not keyword_set (uncalibrated)) then begin
xyouts, 0.92, 0.03, "Calibrated", alignment = 0, /normal
endif else begin
xyouts, 0.92, 0.03, "Uncalibrated", alignment = 0, /normal
endelse
xyouts, 0.91, 0.01, "Instrument", alignment = 1, /normal
xyouts, 0.92, 0.01, platform, alignment = 0, /normal
; print out version information
xyouts, 0.14, 0.01, "IDL Version:", alignment = 1, /normal
xyouts, 0.15, 0.01, idlVersion, alignment = 0, /normal
xyouts, 0.14, 0.03, "Plot Data Version:", alignment = 1, /normal
xyouts, 0.15, 0.03, finVersion + ', ' + finDate, alignment = 0, /normal
xyouts, 0.14, 0.05, "Langley Data Version:", alignment = 1, /normal
xyouts, 0.15, 0.05, flangVersion + ', ' + flangDate, $
alignment = 0, /normal
xyouts, 0.14, 0.07, "Quicklook Version:", alignment = 1, /normal
xyouts, 0.15, 0.07, qlVersion, alignment = 0, /normal
; write out png file
if (!D.name eq 'Z') then begin
get_quicklook_dir, dateYYMMDD, uncal + platform + 'langley', site, $
facility, 'c1', /create, $
quicklook_dir = ql_dir, errornum = errnum
ql_dir = ql_dir + '/'
fname = ql_dir + anal_ql + ".ch" + $
string(format='(I1.1)', channel) + '.' + $
Tstamp + ".png"
if (n_elements(quiet) eq 0) then printf, zlog, $
"Output file is " + fname
tvlct, r, g, b, /get
write_png, fname, tvrd(), r, g, b
endif
;; Close this file, if we are going to open a new one for the
;; next plot
if (!D.name eq 'PS' and n_elements(onefile) eq 0) then device, /close
;; If we are plotting to screen, this lets us prompt to move on
foo=''
if (!D.name eq 'X') then read,foo,prompt='Hit return to continue'
;; used if /noplot is set - see above.
skip_plot:
endfor ;; channel loop
endfor ;; file loop
if (!D.name eq 'PS') then device, /close
;; print out our statistics, if /stats was set
if (n_elements(stats) gt 0) then begin
if (td ne '') then printf, zlog, td, ":"
printf, zlog, "Channel ", [0,1,2,3,4,5,6]
printf, zlog, "# of plots", nplots
printf, zlog, "Accepted ", badflag_counter(*,0)
printf, zlog, "Badflag=1 ", badflag_counter(*,1)
printf, zlog, "Badflag=2 ", badflag_counter(*,2)
printf, zlog, "Badflag=3 ", badflag_counter(*,3)
endif
; close the log file and print to the zebra logger
if (zlog eq 2) then begin
printf, zlog, "Langley quicklook finished normally"
printf, zlog, " "
close, zlog
spawn, zlog_command
endif
; reset the output to the starting output
set_plot, Dname
end
src/langley_batch.pro 0 → 100644
View file @ 02287497
; $Id: langley_batch.pro,v 1.9 2012-01-09 16:03:19 koontz Exp $
;--------------------------------------------------------------------------------
;+
; Abstract:
;
; Author:
; Todd Halter, PNNL
;
; Date Created:
; November, 1999
;
; Date Last Modified:
; $Date: 2012-01-09 16:03:19 $
;
;-
pro langley_batch
date = long (getenv ("DATE"))
platform = getenv ("PLATFORM")
site = getenv ("SITE")
facility = getenv ("FACILITY")
algorithm = getenv ("ALGORITHM")
resolve_routine, 'langley'
resolve_all
langley, 1, to=6, date = date, platform = platform, site=site, $
facility=facility, algorithm = 'barnard', $
plot_to = 'png'
langley, 1, to=6, date = date, platform = platform, site=site, $
facility=facility, algorithm = 'michalsky', $
plot_to = 'png'
; langley, 1, to=6, date = date, platform = platform, site=site, $
; facility=facility, algorithm = 'barnard', /onefile, $
; plot_to = 'ps', /zlog
; langley, 1, to=6, date = date, platform = platform, site=site, $
; facility=facility, algorithm = 'michalsky', /onefile, $
; plot_to = 'ps', /zlog
return
end
src/langley_retriever.h 0 → 100755
View file @ 02287497
#ifndef _LANGLEY_RETRIEVER_H
#define _LANGLEY_RETRIEVER_H
#include "bw_adi.h"
#ifndef NO_RETRIEVER_H
const char *Platform_List[]=
{
//"nimfr.b1|mfrsr.b1",
"AUTO",
NULL
};
const char *Field_List[][BW_MAX_FIELDS]=
{
{
//mfrsr.b1
"direct_normal_broadband",
"direct_normal_narrowband_filter1",
"direct_normal_narrowband_filter2",
"direct_normal_narrowband_filter3",
"direct_normal_narrowband_filter4",
"direct_normal_narrowband_filter5",
"direct_normal_narrowband_filter6",
"qc_direct_normal_broadband",
"qc_direct_normal_narrowband_filter1",
"qc_direct_normal_narrowband_filter2",
"qc_direct_normal_narrowband_filter3",
"qc_direct_normal_narrowband_filter4",
"qc_direct_normal_narrowband_filter5",
"qc_direct_normal_narrowband_filter6",
"nominal_calibration_factor_broadband",
"nominal_calibration_factor_filter1",
"nominal_calibration_factor_filter2",
"nominal_calibration_factor_filter3",
"nominal_calibration_factor_filter4",
"nominal_calibration_factor_filter5",
"nominal_calibration_factor_filter6",
NULL
},
};
#endif // NO_RETRIEVER_H
#endif // _LANGLEY_RETRIEVER_H
src/lfit.c 0 → 100644
View file @ 02287497
/*
* $Id: lfit.c,v 1.2 1997/07/22 16:27:45 d3a230 vap-langley-2.25-0.sol5_10 $
*
* lfit() does a weighted linear fit, returning all sorts of error
* numbers. You can unweight by using NULL for the w array on input.
*
*
* $Log: lfit.c,v $
* Revision 1.2 1997/07/22 16:27:45 d3a230
* Added some sanity checking when we have two or fewer non-zero-weighted
* points. Pathological, of course, but if we don't do it we get math errors
* all the time.
*
* Revision 1.1 1997/07/18 15:58:43 d3a230
* Initial revision
*
*/
static char *rcsid="$Id: lfit.c,v 1.2 1997/07/22 16:27:45 d3a230 vap-langley-2.25-0.sol5_10 $";
static char *rcsstate="$State: vap-langley-2.25-0.sol5_10 $";
#include <math.h>
/*
* Do a linear fit y = a+bx, weighted by w
*
* Right now the siga and sigb assume no measurement uncertainty - i.e.
* we estimate via sampling theory - sigy =~ sqrt(npr*syx2)
*/
void lfit(float *x, float *y, float *w, int n, float *a, float *b,
float *siga, float *sigb, float *sigfit, int *count)
{
int i;
double sumx=0, sumy=0, sumxy=0, sumx2=0, sumy2=0, sumw=0;
double denom, slope, intercept, npr, m, syx2;
m=0;
/* This lets us set w=NULL on input to do an unweighted fit */
#define W (w ? w[i] : 1)
for (i=0; i<n; i++) {
sumx += W*x[i];
sumx2 += W*x[i]*x[i];
sumy += W*y[i];
sumy2 += W*y[i]*y[i];
sumxy += W*x[i]*y[i];
sumw += W;
/* This is a little goofy - m holds the number of non-zero */
/* weighted points, which we need to get our n/(n-2) factor right */
if (w[i] > 0.0) m += 1.0;
}
#undef W
/* simple formulas */
denom = sumw*sumx2 - sumx*sumx;
slope = (sumw*sumxy - sumx*sumy)/denom;
intercept = (sumy*sumx2 - sumx*sumxy)/denom;
/* I'm not sure this is right, especially with the weights */
/* syx2 is standard error of estimate */
syx2 = (sumy2 - intercept*sumy - slope*sumxy)/sumw;
// Actually, I think that is analytically correct - but there sometimes
// gives a roundoff that leads to negative values. I will look for both
// negative and small absolute values later on
/* sampling error */
npr = m > 2 ? m/(m-2) : 1;
*a = intercept;
*b = slope;
// A large negative value is bad, so blow stuff up
if (syx2 < -1e-6) {
if (siga) *siga = 1e8;
if (sigb) *sigb = 1e8;
if (sigfit) *sigfit = 1e8;
if (count) *count = m;
return;
}
// Now, trap for roundoffs - set the min error to 1e-6. This will
// capture really small errors of both signs
if (syx2 < 1e-6) {
syx2 = 1e-6;
}
// Now calculate errors of the fit
if (siga) *siga = m>2 ? sqrt(syx2*npr*sumx2/denom) : 1e8;
if (sigb) *sigb = m>2 ? sqrt(syx2*npr*sumw/denom) : 1e8;
if (sigfit) *sigfit = m>2 ? sqrt(syx2) : 1e8;
if (count) *count = m;
return;
}
char *get_fit_id()
{
return(rcsid);
}
char *get_fit_state()
{
return(rcsstate);
}
src/michalsky_analysis.c 0 → 100644
View file @ 02287497
#define MICHALSKY_ANALYSIS_C
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <memory.h>
#include <string.h> /* memset here in djgpp */
#include "michalsky_langley.h" /* added by tdh for pnnl code */
#include "la.h"
#include "sunae.h"
/*
** Do langley analysis for a set of observations. The input records should
** be partitioned into whatever conditions are required.
** 'cfg' points to a description of the input data
** 'obs' points to a list of cfg->n_obs input records
**
** the return value is a pointer to a list of optical depth/i0 pairs; the
** number of pairs in the list is determined by how many 'columns' have
** langley analysis requested in cfg->ana.
**
** 6/27/97 rs
** Also count the number of records that are within the specified airmass
** range and return it within the la_return structure.
**
** 7/1/97 rs
** Add first_am and last_am to the items in the la_return structure; these
** are the first and last record times that are in the specified airmass
** range.
**
** 9/4/97 rs
** The data structure that gets returned now includes the points that
** make up the langley event as well as the flags indicating why or
** why not they were used.
**
*/
/* tdh - moved to michalsky_langley.h sturct solar_t */
/* tdh - moved to michalsky_langley.h solar_t *solar_geom */
/* tdh - moved to michalsky_langley.h double soldist2; */
typedef enum
{
Good=0,
LowSunAngle,
OutsideAM,
BeamTooLow,
Delta1,
Delta2,
LSPass1,
LSPass2,
TooFewPts,
SDTooBig
} rec_stat_t;
/* tdh - moved to michalsky_langley.h char **rec_status; */
typedef struct {
double x, y;
int idx;
} pt_idx_t;
pt_idx_t *get_good_pts(la_cfg *, la_obs *, int, int *);
int outside_am(la_cfg *, double);
extern float michalsky_sdist;
la_return *lang_ana(la_cfg *cfg, la_obs *obs)
{
void print_ana(la_cfg *, la_obs *),
init_status(int, int),
check_airmass(la_cfg *, int),
check_magnitude(la_cfg *, la_obs *, int),
remove_cloud_pass(la_cfg *, la_obs *, int),
check_concavity(la_cfg *, la_obs *, int),
remove_sd_outliers(la_cfg *, la_obs *, int),
final_checks(la_cfg *, la_obs *, int),
copy_key_flag(la_cfg *, int),
find_first_last_recs(la_cfg *, la_obs *, la_return *),
copy_points(la_cfg *, la_obs *, la_return *),
count_am(la_cfg *, la_return *),
free_status(la_cfg *);
la_return *lrv,
*lin_regress(la_cfg *, la_obs *);
int chan, n, lo, hi,
key_chan(la_cfg *);
void get_solar_geom(la_cfg *, la_obs *);
init_status(cfg->n_obs, cfg->n_data);
get_solar_geom(cfg,
obs);
if ((n = key_chan(cfg)) >= 0)
{
lo = n;
hi = n + 1;
}
else
{
lo = 0;
hi = cfg->n_data;
}
for (chan=lo; chan < hi; chan++)
{
check_magnitude(cfg, obs, chan);
check_airmass(cfg, chan);
remove_cloud_pass(cfg, obs, chan);
check_concavity(cfg, obs, chan);
remove_sd_outliers(cfg, obs, chan);
final_checks(cfg, obs, chan);
}
if (n >= 0)
{
copy_key_flag(cfg, n);
}
/* lin_regress will figure out the key channel for itself */
lrv = lin_regress(cfg, obs);
count_am(cfg, lrv);
find_first_last_recs(cfg, obs, lrv);
copy_points(cfg, obs, lrv);
print_ana(cfg, obs);
/* tdh - add this function to be able to get the data for a netcdf file */
get_michalsky_data (cfg, obs, lrv);
free(solar_geom);
free_status(cfg);
return lrv;
}
void copy_key_flag(la_cfg *cfg, int n)
{
int i, j, fl;
j = cfg->n_data;
for (i=0; i<cfg->n_obs; i++) {
fl = rec_status[i][n];
memset(rec_status[i], fl, j);
}
}
void free_status(la_cfg *cfg)
{
int i;
for (i=0; i<cfg->n_obs; i++)
free(rec_status[i]);
free(rec_status);
}
void init_status(int length, int width)
{
int i;
if ((rec_status = malloc(length*sizeof(char *))) == NULL) {
fprintf(stderr, "init_status: alloc fail 1 (%d)\n", length);
exit(1);
}
for (i=0; i<length; i++) {
if ((rec_status[i] = malloc(width)) == NULL) {
fprintf(stderr, "init_status: alloc fail 2 (%d)\n", width);
exit(1);
}
memset(rec_status[i], '\0', width);
}
}
void get_solar_geom(la_cfg *cfg,
la_obs *obs)
{
ae_pack aep;
int i, j;
float sdist_avg;
int sdist_count;
if ((solar_geom = malloc(cfg->n_obs*sizeof(solar_t))) == NULL)
{
fprintf(stderr, "get_solar_geom: alloc fail (%ld)\n",
cfg->n_obs*sizeof(solar_t));
exit(1);
}
aep.lat = cfg->lat;
aep.lon = cfg->lon;
aep.is_tst = 0;
sdist_avg = 0.0;
sdist_count = 0;
for (i=0; i<cfg->n_obs; i++)
{
aep.year = obs[i].obs_time.year;
aep.doy = obs[i].obs_time.doy;
aep.hour = obs[i].obs_time.gmt;
sunae(&aep);
solar_geom[i].zen_r = M_DTOR*(90.0 - aep.el);
solar_geom[i].az_r = M_DTOR*aep.az;
sdist_avg = sdist_avg + aep.soldst;
sdist_count++;
soldist2 = aep.soldst*aep.soldst;
if (aep.el > 1.0)
{
solar_geom[i].am = armass(aep.el);
}
else
{
solar_geom[i].am = 0.0;
for (j=0; j<cfg->n_data; j++)
{
rec_status[i][j] = LowSunAngle;
}
}
}
if (sdist_count > 0)
{
michalsky_sdist = sdist_avg / sdist_count;
if (isnan(michalsky_sdist))
{
michalsky_sdist = -9999.;
}
}
else
{
michalsky_sdist = -9999.;
}
}
int key_chan(la_cfg *cfg)
{
int i;
for (i=0; i<cfg->n_data; i++)
if (cfg->ana[i] == Key_Ana)
return i;
return -1;
}
void find_first_last_recs(la_cfg *cfg, la_obs *obs, la_return *l)
{
int i;
for (i=0; i<cfg->n_data; i++) {
l[i].first_am = 0;
l[i].last_am = 0;
}
for (i=0; i<cfg->n_obs && outside_am(cfg, solar_geom[i].am); i++) {
}
if (i == cfg->n_obs)
return;
l->first_am = obs[i].obs_time.tt;
l->last_am = obs[i + l->n_am - 1].obs_time.tt;
for (i=1; i<cfg->n_data; i++) {
l[i].first_am = l[0].first_am;
l[i].last_am = l[0].last_am;
}
}
void check_airmass(la_cfg *cfg, int chan)
{
int i;
for (i=0; i<cfg->n_obs; i++)
if (rec_status[i][chan] == Good)
if (outside_am(cfg, solar_geom[i].am))
rec_status[i][chan] = OutsideAM;
}
void count_am(la_cfg *cfg, la_return *l)
{
int i, n;
for (n=0, i=0; i<cfg->n_obs; i++)
if (!outside_am(cfg, solar_geom[i].am))
n++;
for (i=0; i<cfg->n_data; i++)
l[i].n_am = n;
}
void check_magnitude(la_cfg *cfg, la_obs *obs, int chan)
{
int i;
for (i=0; i<cfg->n_obs; i++)
if (rec_status[i][chan] == Good)
if (obs[i].data[chan] < 0.02)
rec_status[i][chan] = BeamTooLow;
}
void print_ana(la_cfg *cfg, la_obs *obs)
{
int i, j;
extern arg_t arguments;
if (arguments.debug == 0)
return;
if (arguments.debug > 1) {
printf("lat=%.3f\nlon=%.3f\n", cfg->lat, cfg->lon);
printf("n_data=%d\nn_obs=%d\n", cfg->n_data, cfg->n_obs);
printf("Ana:");
for (i=0; i<cfg->n_data; i++)
printf(" %d", cfg->ana[i]);
puts("");
}
for (i=0; i<cfg->n_obs; i++) {
printf("%d %d %.6f", obs[i].obs_time.year, obs[i].obs_time.doy,
obs[i].obs_time.gmt);
printf(" %.2f %.2f %.4f", M_RTOD*solar_geom[i].zen_r,
M_RTOD*solar_geom[i].az_r, solar_geom[i].am);
printf(" ");
for (j=0; j<cfg->n_data; j++)
printf("%d", rec_status[i][j]);
for (j=0; j<cfg->n_data; j++)
printf(" %.3f", obs[i].data[j]);
puts("");
}
}
la_return *lin_regress(la_cfg *cfg, la_obs *obs)
{
la_return *rv, lsf(pt_idx_t *, int);
pt_idx_t *pts;
int i, j, n, key;
if ((rv = malloc(cfg->n_data*sizeof(la_return))) == NULL)
{
fprintf(stderr, "lin_regress: alloc fail 1\n");
return NULL;
}
for (j=0; j<cfg->n_data; j++)
{
if ((pts = get_good_pts(cfg, obs, j, &n)) == NULL)
{
rv[j].op_depth = -9999.;
rv[j].i0 = -9999.;
rv[j].i0_sdist = -9999.;
rv[j].n = 0;
rv[j].sd = -9999.;
continue;
}
rv[j] = lsf(pts, n);
free(pts);
}
/*
** now do the last check - magnitude of standard deviation around
** the least squares line. if there is a key channel, then use
** that channel to set all the cahnnels; otherwise, check each
** individually
*/
if ((key = key_chan(cfg)) >= 0)
{
if (rv[key].sd > cfg->lsfit_sd)
{
j = cfg->n_data;
for (i=0; i<cfg->n_obs; i++)
{
if (rec_status[i][key] == Good)
{
memset(rec_status[i], SDTooBig, j);
}
}
for (i=0; i<cfg->n_data; i++)
{
rv[i].n = 0;
}
}
}
else
{
for (j=0; j<cfg->n_data; j++)
{
if (rv[j].sd > cfg->lsfit_sd)
{
for (i=0; i<cfg->n_obs; i++)
{
if (rec_status[i][j] == Good)
{
rec_status[i][j] = SDTooBig;
}
}
rv[j].n = 0;
}
}
}
return rv;
}
la_return lsf(pt_idx_t *pt, int n)
{
double t1, t2, m, mean_x, mean_y, sxx, sxy, syy;
int i;
la_return rv;
extern arg_t arguments;
if (n < 4)
{
rv.op_depth = 0.0;
rv.i0 = 0.0;
rv.i0_sdist = 0.0;
rv.n = n;
if (arguments.debug > 1)
printf("lsf: too few points to fit (%d)\n", n);
return rv;
}
mean_x = 0.0;
mean_y = 0.0;
for (i=0; i<n; i++)
{
mean_x += pt[i].x;
mean_y += pt[i].y;
}
mean_x /= n;
mean_y /= n;
sxx = 0.0;
sxy = 0.0;
syy = 0.0;
for (i=0; i<n; i++)
{
t1 = pt[i].x - mean_x;
t2 = pt[i].y - mean_y;
sxy += t1*t2;
sxx += t1*t1;
syy += t2*t2;
}
if (isnan (sxy))
{
/* Can't compute if sxy is a NaN */
rv.op_depth = -9999.;
m = -9999.;
rv.i0 = -9999.; /* intercept */
rv.i0_sdist = -9999.; /* correct for solar dist */
rv.sd = -9999.; /* std deviation */
rv.n = -9999.;
}
else
{
if (sxx != 0.0)
{
rv.op_depth = sxy/sxx;
m = rv.op_depth; /* slope */
rv.i0 = exp(mean_y - rv.op_depth*mean_x); /* intercept */
rv.i0_sdist = rv.i0*soldist2; /* correct for solar dist */
rv.sd = sqrt((syy - m*m*sxx)/(n - 2)); /* std deviation */
rv.n = n;
// trs 8/11/15: switch on valid i0, using the valid range of
// 0.5,5000, to prevent outputing infinities or other garbage
if (rv.i0 > 5000. || rv.i0 < 0.05 ) {
rv.op_depth = -9999.;
m = -9999.;
rv.i0 = -9999.; /* intercept */
rv.i0_sdist = -9999.; /* correct for solar dist */
rv.sd = -9999.; /* std deviation */
rv.n = -9999.;
}
}
else
{
/* Can't divide by zero, so set results to -9999. */
rv.op_depth = -9999.;
m = -9999.;
rv.i0 = -9999.; /* intercept */
rv.i0_sdist = -9999.; /* correct for solar dist */
rv.sd = -9999.; /* std deviation */
rv.n = -9999.;
}
}
if (arguments.debug > 1) {
printf("y = %.4fx + %.4f on %d points\n", rv.op_depth, rv.i0, n);
for (i=0; i<n; i++)
printf(" %.4f %.4f\n", pt[i].x, pt[i].y);
}
return rv;
}
/*
** run through the observations twice, and on each pass...
** - compute the least squares fit
** - remove points that are > 1.5x standard deviation from ls line
*/
void remove_sd_outliers(la_cfg *cfg, la_obs *obs, int chan)
{
la_return rv, lsf(pt_idx_t *, int);
pt_idx_t *pts;
int pass, n;
void check_sd(la_cfg *, la_obs *, la_return *, int, int);
for (pass=0; pass<2; pass++)
{
if ((pts = get_good_pts(cfg, obs, chan, &n)) == NULL)
{
continue;
}
rv = lsf(pts, n);
check_sd(cfg, obs, &rv, chan, pass);
free(pts);
}
}
/*
** loop through the observations; compute the distance from the
** least squares line for channel 'j' in 'rv'; if more than 1.5*sd,
** then mark the observation as bad1 or bad2, depending on 'pass'
*/
void check_sd(la_cfg *cfg, la_obs *obs, la_return *rv, int j, int pass)
{
int i;
double sd, y, dif, i0;
sd = cfg->out_limit*rv->sd;
/* js - added a check for bad rv->io values. */
if(rv->i0 <= 0) {
for (i=0; i<cfg->n_obs; i++) {
if (rec_status[i][j] == Good)
rec_status[i][j] = pass == 0 ? LSPass1 : LSPass2;
}
return;
}
i0 = log(rv->i0);
for (i=0; i<cfg->n_obs; i++) {
if (rec_status[i][j] != Good)
continue;
/* the lsf value for the current obs */
y = rv->op_depth*solar_geom[i].am + i0;
dif = fabs(y - log(obs[i].data[j]));
if (dif > sd)
rec_status[i][j] = pass == 0 ? LSPass1 : LSPass2;
}
}
/*
** this checks for and tries to remove periods of time when clouds
** block the sun; the scheme is to check for positive slope in the plot
** of airmass vs. intensity; this should be physically impossible if
** the sky is really clear, and so indicates "clearing". once you've
** decided the the sun is getting brighter, go back and find the
** darkest recent record, then remove around it.
**
*/
void remove_cloud_pass(la_cfg *cfg, la_obs *obs, int chan)
{
void rm_clouds(pt_idx_t *, int, la_cfg *, la_obs *, int);
pt_idx_t *pts;
int n,
pt_cmp(const void *, const void *);
if ((pts = get_good_pts(cfg, obs, chan, &n)) == NULL)
return;
qsort(pts, n, sizeof(pt_idx_t), pt_cmp);
rm_clouds(pts, n, cfg, obs, chan);
free(pts);
}
/*
** 'pts' is sorted by airmass; if intensity increases, then a cloud
** is passing; make a mark (m1) at the first point of increase, then
** keep moving forward until intensity decreases again (m2); then
** mark all records between time m1 +/- m2 bad.
*/
void rm_clouds(pt_idx_t *pts, int n, la_cfg *cfg, la_obs *obs, int chan)
{
int i, m1, m2;
time_t t1, t2;
void clobber(la_cfg *, la_obs *, time_t, time_t, int, rec_stat_t);
extern arg_t arguments;
m1 = m2 = -1;
for (i=0; i<n-1; i++) {
if (arguments.debug > 1)
printf("%d %.5f %.5f %d %d\n", i, pts[i].x, pts[i].y, m1, m2);
/* not sure if diff_slop is being applied properly -- shouldn't it be +- diff_slop ? */
if (pts[i].y + cfg->diff_slop < pts[i+1].y) {
if (m1 < 0)
m1 = i;
} else {
if (m1 >= 0) {
m2 = i;
if (arguments.debug > 1)
printf("%d %d clobbered\n", m1, m2);
t2 = obs[pts[m2].idx].obs_time.tt;
t1 = obs[pts[m1].idx].obs_time.tt;
clobber(cfg, obs, t1, t2, chan, Delta1);
m1 = -1;
}
}
}
if (m1 != -1) { /* positive slope on last 2 points */
t1 = obs[pts[m1].idx].obs_time.tt;
t2 = obs[pts[n-1].idx].obs_time.tt;
clobber(cfg, obs, t1, t2, chan, Delta1);
}
}
void clobber(la_cfg *cfg, la_obs *obs, time_t base,
time_t end, int chan, rec_stat_t st)
{
int i;
time_t t1, t2;
if (end > base) {
t2 = end;
t1 = base - (end - base);
} else {
t1 = end;
t2 = base + (base - end);
}
/*
** the next bit of trickery is kind of arbitrary, but is intended
** to add a bit of slop to the cloberred area...the rationale is
** that the ASCII representation of the time of day may be imprecise,
** so add a bit more to the beginning and end of the "hole."
*/
t1 -= 3; /* seconds */
t2 += 3;
for (i=0; i<cfg->n_obs; i++)
if (rec_status[i][chan] == Good &&
obs[i].obs_time.tt >= t1 && obs[i].obs_time.tt <= t2) {
rec_status[i][chan] = st;
}
/*
for (i=0; i<cfg->n_obs; i++)
if (obs[i].obs_time.tt >= t1 && obs[i].obs_time.tt <= t2)
if (rec_status[i][chan] == Good)
rec_status[i][chan] = st;
*/
}
int pt_cmp(const void *v1, const void *v2)
{
pt_idx_t *p1 = v1, *p2 = v2;
if (p1->x < p2->x)
return -1;
else if (p1->x > p2->x)
return 1;
else
return 0;
}
/*
** this filter is intended to check the second derivative; it's implementation
** is maybe not exactly the same as Lee's.
*/
void check_concavity(la_cfg *cfg, la_obs *obs, int chan)
{
int n,
pt_cmp(const void *, const void *);
pt_idx_t *pts;
void deriv2(la_cfg *, la_obs *, pt_idx_t *, int, int);
if ((pts = get_good_pts(cfg, obs, chan, &n)) == NULL)
return;
qsort(pts, n, sizeof(pt_idx_t), pt_cmp);
deriv2(cfg, obs, pts, n, chan);
free(pts);
}
/*
** loop through each pair of points in the list; compute slope
** between each; compute the average of all of them;
** Go through the list again; if slope between 2 points is great
** than 2 times the average slope, throw away the 2nd of the 2.
**
** I think.
*/
void deriv2(la_cfg *cfg, la_obs *obs, pt_idx_t *pts, int n, int chan)
{
double *slopes, avg;
int i;
if (n < 3)
return;
if ((slopes = malloc((n-1)*sizeof(double))) == NULL) {
fprintf(stderr, "deriv2: alloc fail %d doubles\n", n-1);
exit(1);
}
avg = 0.0;
for (i=0; i<n-1; i++) {
slopes[i] = (pts[i+1].y - pts[i].y) / (pts[i+1].x - pts[i].x);
avg += slopes[i];
}
avg /= n - 1;
#ifdef MEANDERIV
/*
** compare each slope against the average of all the slopes
*/
for (i=0; i<n-1; i++)
if (slopes[i] < 2*avg)
rec_status[pts[i+1].idx][chan] = Delta2;
#else
/*
** compare each slope to the previous
*/
for (i=0; i<n-2; i++)
if (slopes[i] < 2*slopes[i+1])
rec_status[pts[i+1].idx][chan] = Delta2;
#endif
free(slopes);
}
/*
** make sure that at least 1/3 of the original number of points still remain
** - compare to the number in 2..6 airmasses...
*/
void final_checks(la_cfg *cfg, la_obs *obs, int chan)
{
int i, n, n_in;
n_in = 0;
for (i=0; i<cfg->n_obs; i++)
if (solar_geom[i].am >= cfg->low_am && solar_geom[i].am <= cfg->high_am)
n_in++;
n = 0;
for (i=0; i<cfg->n_obs; i++)
if (rec_status[i][chan] == Good)
n++;
if (n < cfg->frac_pts*n_in)
for (i=0; i<cfg->n_obs; i++)
if (rec_status[i][chan] == Good)
rec_status[i][chan] = TooFewPts;
}
/*
** return a list of the observations for channel 'chan' which are still
** marked 'good'.
*/
pt_idx_t *get_good_pts(la_cfg *cfg, la_obs *obs, int chan, int *got)
{
pt_idx_t *pts;
int i, n;
for (n=0, i=0; i<cfg->n_obs; i++)
if (rec_status[i][chan] == Good)
n++;
if (n == 0) {
*got = 0;
return NULL;
}
if ((pts = malloc(n*sizeof(pt_idx_t))) == NULL) {
fprintf(stderr, "get_good_pts: alloc fail %d\n", n);
exit(1);
}
for (n=0, i=0; i<cfg->n_obs; i++)
if (rec_status[i][chan] == Good) {
pts[n].x = solar_geom[i].am;
pts[n].y = log(obs[i].data[chan]);
pts[n].idx = i;
n++;
}
*got = n;
return pts;
}
/*
** Make lists of the points within the given airmass range; copy the
** x, y, and flag values.
*/
void copy_points(la_cfg *cfg, la_obs *obs, la_return *lrv)
{
int i, j, k;
for (i=0; i<cfg->n_data; i++) {
if (lrv[i].n_am == 0) {
lrv[i].x = NULL;
lrv[i].y = NULL;
continue;
}
if ((lrv[i].x = malloc(lrv[i].n_am * sizeof(double))) == NULL) {
fprintf(stderr, "copy_points: malloc fail (1)\n");
exit(1);
}
if ((lrv[i].y = malloc(lrv[i].n_am * sizeof(pt_t))) == NULL) {
fprintf(stderr, "copy_points: malloc fail (2)\n");
exit(1);
}
j = 0;
for (k=0; k<cfg->n_obs; k++) {
if (outside_am(cfg, solar_geom[k].am))
continue;
lrv[i].x[j] = solar_geom[k].am;
lrv[i].y[j].y = obs[k].data[i];
lrv[i].y[j].flag = rec_status[k][i];
j++;
}
if (j != lrv[i].n_am) {
fprintf(stderr, "Fatal: copy_points: j <> n_am\n");
exit(1);
}
}
}
int outside_am(la_cfg *cfg, double val)
{
return val < cfg->low_am || val > cfg->high_am;
}
#undef MICHALSKY_ANALYSIS_C
src/michalsky_langley.c 0 → 100644
View file @ 02287497
/*******************************************************************************
* COPYRIGHT (C) 2001 Battelle Memorial Institute.
* All Rights Reserved. (Now, with this out of the way...)
*
* RCS INFORMATION:
* $RCSfile: michalsky_langley.c,v $
* $Revision: 1.14 $
* $Author: koontz $
* $Locker: $
* $Date: 2011/12/29 18:57:13 $
* $State: Exp $
* $Name: $
* $Id: michalsky_langley.c,v 1.14 2011/12/29 18:57:13 koontz Exp $
*
* FUNCTIONS IN THIS FILE:
*
* DESIGN:
* <Discuss overall design/purpose of class.>
*******************************************************************************/
#define MICHALSKY_LANGLEY_C
/****** General Includes ******/
#include <stdio.h>
#include <time.h>
#include <math.h>
/****** Zebra Includes ******/
/* #include "defs.h"
* #include "message.h"
* #include "DataStore.h"
*/
/****** BW Includes ******/
/* #define BW_CODE
* #include "bw_main.h"
*/
/****** Application Includes ******/
#include "langley.h"
#include "michalsky_langley.h"
#define BW_CODE
#include "bw_adi.h"
#define TRUE 1
#define FALSE 0
float michalsky_sdist; /* To capture sdist from gest_solar_geom */
/*******************************************************************************
* Function:
* michalsky_langley -
*
* Inputs:
*******************************************************************************/
int michalsky_langley (DATA *D, DATA *newD, float *gNomCal)
{
/* functions used by michalsky_langley (... */
int setup_for_michalsky (DATA *);
int write_newD (DATA *, float *);
int status;
setup_for_michalsky (D);
/* I don't think this is needed
config_file ();*/ /* This is michalsky code */
langley (); /* This is michalsky code */
if ((status = write_newD (newD, gNomCal)) < 0) {
bw_return(-1, "Michalsky: problem writing newD");
}
return (TRUE);
} /* michalsky_langley (... */
int setup_for_michalsky (DATA *D)
{
char cfg_file[512];
int i, o, t, ot, r;
int sot, eot, valid_idx;
rec_t *rec;
extern long startime, etime;
/* functions used by setup_for_michalsky (... */
/*
* init global vars
*/
michalsky_npart = 0;
/* determin the number of samples */
ot = eot = 0;
sot = -1;
for (o=0; o<D->nObs[B1]; o++)
{
for (t=0; t<D->nSamples[B1][o][0][IN_BROADBAND]; t++, ot++)
{
if (D->sampleTimes[B1][o][0][IN_BROADBAND][t].zt_Sec >= startime &&
D->sampleTimes[B1][o][0][IN_BROADBAND][t].zt_Sec <= etime)
{
if (sot == -1) sot = ot;
eot = ot;
} /* end if (D->sampleTimes... */
} /* end for (t... */
} /* end for (o... */
if (sot == -1)
bw_return (-1, "Michalsky: no samples for given time period");
la_n_records = (eot - sot) + 10; /* was +1, changed to +10 because code bombs later ... A. Koontz, PNNL, June, 2013 */
printf ("DEBUG: la_n_records = %d\n", la_n_records);
/* set up arrays to hold michalsky data */
la_records = CALLOC (la_n_records, rec_t);
michalsky_status = CALLOC (la_n_records, int *);
michalsky_time = CALLOC (la_n_records, long);
michalsky_am = CALLOC (la_n_records, float);
michalsky_lnI = CALLOC (la_n_records, float *);
for (r=0; r<la_n_records; r++)
{
michalsky_lnI[r] = CALLOC (NCHANNELS, float);
michalsky_status[r] = CALLOC (NCHANNELS, int);
} /* for (r... */
/*
* set up arrays to hold michalsky return data.
* Using NSLABS as a guess, way to many but will work
* should have 2 obs for every day.
*/
michalsky_ret = CALLOC (NSLABS, la_return *);
for (o=0; o < NSLABS; o++)
michalsky_ret[o] = CALLOC (NCHANNELS, la_return);
/*
* populate data vars
*/
ot = 0;
valid_idx = 0;
for (o=0; o < D->nObs[B1]; o++)
{
for (t=0; t < D->nSamples[B1][o][0][IN_BROADBAND]; t++, ot++)
{
if (ot >= sot && ot <= eot)
{
#ifdef DEBUGIT
printf ("DEBUG: michalsky_langley at AAA: valid_idx = %d\n", valid_idx);
#endif
rec = &la_records[valid_idx];
jtoydm ((time_t) D->sampleTimes[B1][o][0][IN_BROADBAND][t].zt_Sec,
rec);
rec->n_direct = NCHANNELS;
rec->direct = CALLOC (NCHANNELS, double);
for (r=0; r < (NCHANNELS - 1); r++)
{
/* if new file format else old format */
rec->direct[r] = D->BWdata[B1][o][0][IN_CHN1+r][t][0];
}
rec->direct[NCHANNELS-1]= D->BWdata[B1][o][0][IN_BROADBAND][t][0];
valid_idx++;
} /* end if (ot... */
} /* for (t... */
} /* for (o... */
/*
* For the michalsky algorithm, broadband needs to be last
*/
for (i=0; i < NCHANNELS - 1; i++)
channels[i] = channels[i + 1];
channels[NCHANNELS - 1] = 0;
/*
* Setting up the arguments structure
*/
/* tdh sprintf (cfg_file, "%s/langley.cfg", getenv ("VAP_CONF_HOME"));
arguments.cfg_fnm = strdup (cfg_file);
*/
arguments.in_fnm = NULL;
arguments.debug = 0;
arguments.julian_out = 1;
arguments.version = 0;
arguments.smooth = 0;
arguments.parse_output = 1;
arguments.p_type = None;
/*
* set up other needed globals
*/
hours_gmt = 0; /* set time in gmt */
return (TRUE);
} /* setup_for_michalsky (... */
int write_newD (DATA *newD, float *gNomCal)
{
int o, f, t, ot, r, valid_idx;
int nobs[NOUTPLATS], chn;
long midpoint, timediff;
/* functions */
/* newD structure is setup and initialized in the barnard algroithm. */
/*
* Populate LANGLEY
* Only has 1 obs
*/
/*
* if the barnard and michalsky algorithms do not produce the same number
* of data points exit out
*/
if (newD->nSamples[LANGLEY][0][0][OD] != michalsky_npart)
{
bw_return (-1, "The algorithms do not give the same number of langleys.");
}
for (t=0; t<michalsky_npart; t++)
{
/* write out the sample time as the midpoint in the data used */
midpoint = ((michalsky_ret[t][0].last_am - michalsky_ret[t][0].first_am) /
2) + michalsky_ret[t][0].first_am;
if (labs (midpoint - newD->sampleTimes[LANGLEY][0][0][OD][t].zt_Sec) > 3600)
{
bw_return (-1, "Algorithm output times do not match within 1hr.");
}
for (r=0; r < NCHANNELS; r++)
{
/*
* the barnard chn flds come first then the then michalsky
* barnard fld start + NCHANNELS = michalsky fld start
*/
/* broadband need to be first in output */
chn = (r == 0) ? NCHANNELS - 1 : r - 1;
if ((int) michalsky_ret[t][chn].op_depth == -9999)
{
newD->BWdata[LANGLEY][0][0][OD + NCHANNELS + r][t][0] = -9999.;
}
else
{
newD->BWdata[LANGLEY][0][0][OD + NCHANNELS + r][t][0] =
(michalsky_ret[t][chn].op_depth == 0) ?
0 : -michalsky_ret[t][chn].op_depth;
}
/* Convert the Io values back to "counts" (undo nominal calibration done by ingest */
newD->BWdata[LANGLEY][0][0][SC + NCHANNELS + r][t][0] =
michalsky_ret[t][chn].i0;
if (michalsky_ret[t][chn].i0 > 0.0)
{
newD->BWdata[LANGLEY][0][0][SC + NCHANNELS + r][t][0] =
newD->BWdata[LANGLEY][0][0][SC + NCHANNELS + r][t][0] * gNomCal[r];
}
/* barnard does not have a SC_NORM */
newD->BWdata[LANGLEY][0][0][SC_NORM + r][t][0] =
michalsky_ret[t][chn].i0_sdist;
if (michalsky_ret[t][chn].i0_sdist > 0.0)
{
newD->BWdata[LANGLEY][0][0][SC_NORM + r][t][0] =
newD->BWdata[LANGLEY][0][0][SC_NORM + r][t][0] * gNomCal[r];
}
/* barnard does not have a SD */
newD->BWdata[LANGLEY][0][0][SD + r][t][0] =
michalsky_ret[t][chn].sd;
newD->BWdata[LANGLEY][0][0][NPTS + NCHANNELS + r][t][0] =
michalsky_ret[t][chn].n;
newD->BWdata[LANGLEY][0][0][GOOD_FRACT + NCHANNELS + r][t][0] =
(float)michalsky_ret[t][chn].n / (float)michalsky_ret[t][chn].n_am;
if (michalsky_ret[t][chn].op_depth == 0.0)
newD->BWdata[LANGLEY][0][0][BAD_FLAG + NCHANNELS + r][t][0] = 1;
else
{
if (michalsky_ret[t][chn].n == 0.0)
newD->BWdata[LANGLEY][0][0][BAD_FLAG + NCHANNELS + r][t][0] = 2;
else
newD->BWdata[LANGLEY][0][0][BAD_FLAG + NCHANNELS + r][t][0] = 0;
} /* end else */
} /* end for (r... */
} /* end for (t... */
/* for michalsky sun to earth distance ... may need to average ??? */
newD->BWdata[LANGLEY][0][0][MSDIST][0][0] = michalsky_sdist;
/*
* Populate PLOT
*/
ot = 0;
for (o=0; o<michalsky_npart; o++)
{
if (michalsky_npart != newD->nObs[PLOT])
{
/* set number of samples = 0 for o */
for (f=0; f < newD->nFields[PLOT][o]; f++)
newD->nSamples[PLOT][o][0][f] = 0;
msg_ELog (EF_INFO, "Plot data did not match for obs %i.", o);
return (FALSE);
} /* end if (michalsky_npart...*/
valid_idx = 0;
for (t=0; t < michalsky_nobs_part[o]; t++)
{
if (michalsky_time[ot] >=
newD->sampleTimes[PLOT][o][0][AM][0].zt_Sec &&
valid_idx <= newD->nSamples[PLOT][o][0][AM] - 1)
{
timediff = newD->sampleTimes[PLOT][o][0][AM][valid_idx].zt_Sec -
michalsky_time[ot];
if (timediff != 0.0)
msg_ELog (EF_INFO, "Sample times do not match for am: %f off by %ds",
michalsky_am[ot], timediff);
for (r=0; r<NCHANNELS; r++)
{
/* broadband need to be first in output */
chn = (r == 0) ? NCHANNELS - 1 : r - 1;
newD->BWdata[PLOT][o][0][LN_I + NCHANNELS + r][valid_idx][0] =
michalsky_lnI[ot][chn];
newD->BWdata[PLOT][o][0][REJECTED + NCHANNELS + r][valid_idx][0] =
michalsky_status[ot][chn];
} /* end for (r... */
newD->BWdata[PLOT][o][0][AM + 1][valid_idx][0] = michalsky_am[ot];
valid_idx++;
} /* end if (michalsky_am[ot]... */
ot++;
} /* end for (t=0... */
if (valid_idx != newD->nSamples[PLOT][o][0][LN_I])
msg_ELog (EF_INFO, "Number of plot samples for obs %i does not match",
o);
} /* end for (o=0... */
return (TRUE);
} /* write_newD (... */
int jtoydm (const time_t ztime,
rec_t *rec)
{
struct tm *ts;
#ifdef DEBUGIT
printf ("DEBUG jtoydm: ztime = %d\n", ztime);
#endif
ts = gmtime (&ztime);
rec->year = ts->tm_year;
rec->doy = ts->tm_yday + 1; /* tm_yday starts at 0 */
rec->tod = ts->tm_hour + ts->tm_min/60.0 + ts->tm_sec/3600.0;
#ifdef DEBUGIT
printf ("DEBUG jtoydm: year=%d doy=%d tod=%lf\n", rec->year, rec->doy, rec->tod);
#endif
return TRUE;
} /* jtoydm (... */
int
get_michalsky_data (la_cfg *cfg, la_obs *obs, la_return *l)
{
static int idx = 0;
static int o = 0;
int t, r;
/*
* if there are not any obs in the valid am range then exit
*/
if (l[0].first_am == 0 || l[0].last_am == 0)
return FALSE;
/*
* Data for langley
*/
for (r=0; r<cfg->n_data; r++)
{
michalsky_ret[o][r].op_depth = l[r].op_depth;
michalsky_ret[o][r].i0 = l[r].i0;
michalsky_ret[o][r].i0_sdist = l[r].i0_sdist;
michalsky_ret[o][r].sd = l[r].sd;
michalsky_ret[o][r].n = l[r].n;
michalsky_ret[o][r].n_am = l[r].n_am;
michalsky_ret[o][r].first_am = l[r].first_am;
michalsky_ret[o][r].last_am = l[r].last_am;
michalsky_ret[o][r].x = NULL;
michalsky_ret[o][r].y = NULL;
} /* end for (r... */
o++;
/*
* Data for langplot
*/
michalsky_nobs_part[michalsky_npart++] = cfg->n_obs;
for (t=0; t<cfg->n_obs; t++, idx++)
{
michalsky_time[idx] = obs[t].obs_time.tt;
michalsky_am[idx] = solar_geom[t].am;
for (r=0; r<cfg->n_data; r++)
{
/* this will pick up bad - values and missing values */
// Rarely, we get obs[t].data pointing to NULL, caused by our ten
// sample buffer on the end (cf. line 116, for la_n_records). So we
// need to make sure our pointer is non-null before using it. -trs, 6/15/15
if (! obs[t].data || obs[t].data[r] <= 0.0)
michalsky_lnI[idx][r] = 0.0;
else
michalsky_lnI[idx][r] = log (obs[t].data[r] * soldist2);
michalsky_status[idx][r] = ((int)rec_status[t][r] == 0) ? 0 : 1;
} /* end for (r... */
} /* end for (t... */
return TRUE;
} /* get_michalsky_data (... */
#undef MICHALSKY_LANGLEY_C
/* michalsky_langley.c */
src/michalsky_langley.h 0 → 100644
View file @ 02287497
/*******************************************************************************
* COPYRIGHT (C) 2001 Battelle Memorial Institute.
* All Rights Reserved. (Now, with this out of the way...)
*
* RCS INFORMATION:
* $RCSfile: michalsky_langley.h,v $
* $Revision: 1.6 $
* $Author: koontz $
* $Locker: $
* $Date: 2011/12/29 18:57:21 $
* $State: Exp $
* $Name: $
* $Id: michalsky_langley.h,v 1.6 2011/12/29 18:57:21 koontz Exp $
*
* LIBRARY DEPENDENCIES:
*
* DESCRIPTION:
* <Discuss overall design/purpose of class.>
*******************************************************************************/
#ifndef MICHALSKY_LANGLEY_H
#define MICHALSKY_LANGLEY_H
/****** system includes ******/
/****** general includes ******/
/****** library includes ******/
/****** application includes ******/
#include "langley.h"
/*
* sets la code to be run as a module
*/
#define WITH_TU
/*
* set to compare each slope in check_concavity(...) against the average
* of all the slopes
* This is set in the michalsky Makefile.
*/
#define MEANDERIV
typedef enum
{
No_Ana, /* don't do langley */
Do_Ana, /* do lanley for corresponding column */
Key_Ana /* this is the Key wavelength */
} ana_t; /* enum ana_t */
enum plot_type
{
None,
Gif,
PS,
Screen
}; /* enum plot_type */
typedef struct
{
int debug;
int parse_output;
int version;
int julian_out;
int smooth;
char *cfg_fnm, in_fnm;
int p_type;
} arg_t;
typedef struct
{
char *var, *val;
int checked;
} field_t;
typedef struct
{
field_t *fields;
int n_fields;
char julian;
} config_t;
typedef struct
{
double lat, lon, /* site lat/lon, degrees */
diff_slop, /* allowable noise in delta filter */
lsfit_sd, /* last filter limit (0.006) */
low_am,
high_am, /* airmass limits */
out_limit, /* how much is an outlier (1.5) */
frac_pts; /* fraction of points after filters */
ana_t *ana; /* list of things to do for each 'column' */
int n_data, /* number of 'columns' per record */
n_obs; /* number of records in group */
} la_cfg;
typedef struct
{
int year, doy, n_direct;
double tod, *direct;
} rec_t;
typedef struct
{
int year, doy;
double gmt;
time_t tt;
} obs_time_t;
typedef struct
{
obs_time_t obs_time; /* time of observation */
double *data; /* list of length n_data data points */
} la_obs;
typedef struct
{
double am, zen_r, az_r;
} solar_t;
typedef struct {
double y;
char flag;
} pt_t;
typedef struct {
double op_depth, /* optical depth */
i0, /* extraterrestrial */
i0_sdist, /* extraterrestrial corrected for solar dist */
sd; /* standard deviation around regression line */
int n, /* number of records used */
n_am; /* number of record in airmass range */
time_t first_am, /* time of first record in am range */
last_am; /* time of last record in am range */
double *x; /* vector of length n_am; airmass for points */
pt_t *y; /* vector of length n_am; */
} la_return;
#endif /* MICHALSKY_LANGLEY_H */
#ifdef MICHALSKY_LANGLEY_C
char **rec_status;
int **michalsky_status, la_n_records, michalsky_npart;
int michalsky_nobs_part[16] = {0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0};
long *michalsky_time;
float **michalsky_lnI, *michalsky_am;
double hours_gmt;
double soldist2;
arg_t arguments;
config_t configuration;
rec_t *la_records = NULL;
solar_t *solar_geom = NULL;
la_return **michalsky_ret = NULL;
#else
extern char **rec_status;
extern int **michalsky_status, la_n_records, michalsky_npart;
extern int michalsky_channels[NCHANNELS];
int michalsky_nobs_part[16];
long *michalsky_time;
extern float **michalsky_lnI, *michalsky_am;
extern double hours_gmt;
extern double soldist2;
extern arg_t arguments;
extern config_t configuration;
extern rec_t *la_records;
extern solar_t *solar_geom;
extern la_return **michalsky_ret;
#endif /* LA_C */
/* function prototypes */
int get_michalsky_data (la_cfg *, la_obs *, la_return *);
la_return *lang_ana(la_cfg *cfg, la_obs *obs);
/* michalsky_langley.h */
src/solarpos.c 0 → 100644
View file @ 02287497
/*
* $Id: solarpos.c 47519 2013-09-10 17:50:19Z shippert $
* Abstract:
* calculates the position of the sun
* Author:
* Nels Larson
* $Log: solarpos.c,v $
* Revision 3.1 2003/03/04 23:49:24 shippert
* Added a smoother to refraction, so it varies linearly from it's value at
* alt == -1 to zero at alt == -2. This eliminates the discontinuity at alt
* == -1, and might even be kinda scientifically justified.
*
* Revision 3.0 1997/07/18 15:55:18 d3a230
* Added rcsid and rcsstate variables and return functions, and set the
* version number to match Nels's (i.e. 3.0).
*
*
* Revision 1.1 1994/06/23 22:37:41 d3a230
* Initial revision
*
* Revision 1.1 1994/05/06 17:40:35 d3a038
* Initial revision
*
*
*/
/* "solarpos.c" contains the integer function solarposition() for calculating
* the position of the Sun as seen from a place on Earth at a specific time.
*
* Version 3.0 - February 20, 1992.
*
* solarposition() employs the low precision formulas for the Sun's coordinates
* given in the "Astronomical Almanac" of 1990 to compute the Sun's apparent
* right ascension, apparent declination, altitude, atmospheric refraction
* correction applicable to the altitude, azimuth, and distance from Earth.
* The "Astronomical Almanac" (A. A.) states a precision of 0.01 degree for the
* apparent coordinates between the years 1950 and 2050, and an accuracy of
* 0.1 arc minute for refraction at altitudes of at least 15 degrees.
*
* The following assumptions and simplifications are made:
* -> refraction is calculated for standard atmosphere pressure and temperature
* at sea level.
* -> diurnal parallax is ignored, resulting in 0 to 9 arc seconds error in
* apparent position.
* -> diurnal aberration is also ignored, resulting in 0 to 0.02 second error
* in right ascension and 0 to 0.3 arc second error in declination.
* -> geodetic site coordinates are used, without correction for polar motion
* (maximum amplitude of 0.3 arc second) and local gravity anomalies.
* -> local mean sidereal time is substituted for local apparent sidereal time
* in computing the local hour angle of the Sun, resulting in an error of
* about 0 to 1 second of time as determined explicitly by the equation of
* the equinoxes.
*
* Right ascension is measured in hours from 0 to 24, and declination in
* degrees from 90 to -90.
* Altitude is measured from 0 degrees at the horizon to 90 at the zenith or
* -90 at the nadir. Azimuth is measured from 0 to 360 degrees starting at
* north and increasing toward the east at 90.
* The refraction correction should be added to the altitude if Earth's
* atmosphere is to be accounted for.
* Solar distance from Earth is in astronomical units, 1 a.u. representing the
* mean value.
*
* The necessary input parameters are:
* -> the date, specified in one of three ways:
* 1) year, month, day.fraction
* 2) year, daynumber.fraction
* 3) days.fraction elapsed since January 0, 1900.
* -> site geodetic (geographic) latitude and longitude.
*
* Refer to the function declaration for the parameter type specifications and
* formats.
*
* solarposition() returns -1 if an input parameter is out of bounds, or 0 if
* values were written to the locations specified by the output parameters.
*/
/* Author: Nels Larson
* Pacific Northwest Lab.
* P.O. Box 999
* Richland, WA 99352
* U.S.A.
*/
static char *rcsid="$Id: solarpos.c 47519 2013-09-10 17:50:19Z shippert $";
static char *rcsstate="$State: process-vap-solarpos-1.1-0 $";
/* Math definitions. */
#define PI 3.1415926535897932
#define TWOPI 6.2831853071795864
#define DEG_RAD 0.017453292519943295
#define RAD_DEG 57.295779513082323
double acos(),
asin(),
atan2(),
cos(),
fabs(),
modf(),
sin(),
tan();
int solarposition(year, month, day, days_1900, latitude, longitude,
ap_ra, ap_dec, altitude, refraction, azimuth, distance)
int year, /* Four digit year (Gregorian calendar).
* [1950 through 2049; 0 o.k. if using days_1900] */
month; /* Month number.
* [1 through 12; 0 o.k. if using daynumber for day] */
double day, /* Calendar day.fraction, or daynumber.fraction.
* [If month is NOT 0:
* 0 through 32; 31st @ 18:10:00 UT = 31.75694
* If month IS 0:
* 0 through 367; 366 @ 18:10:00 UT = 366.75694] */
days_1900, /* Days since 1900 January 0 @ 00:00:00 UT.
* [18262.0 (1950/01/00) through 54788.0 (2049/12/32);
* 1990/01/01 @ 18:10:00 UT = 32873.75694;
* 0.0 o.k. if using {year, month, day} or
* {year, daynumber}] */
latitude, /* Observation site geographic latitude.
* [degrees.fraction, North positive] */
longitude, /* Observation site geographic longitude.
* [degrees.fraction, East positive] */
*ap_ra, /* Apparent solar right ascension.
* [hours; 0.0 <= *ap_ra < 24.0] */
*ap_dec, /* Apparent solar declination.
* [degrees; -90.0 <= *ap_dec <= 90.0] */
*altitude, /* Solar altitude, uncorrected for refraction.
* [degrees; -90.0 <= *altitude <= 90.0] */
*refraction, /* Refraction correction for solar altitude.
* Add this to altitude to compensate for refraction.
* [degrees; 0.0 <= *refraction] */
*azimuth, /* Solar azimuth.
* [degrees; 0.0 <= *azimuth < 360.0, East is 90.0] */
*distance; /* Distance of Sun from Earth (heliocentric-geocentric).
* [astronomical units; 1 a.u. is mean distance] */
{
int daynum(); /* Computes a sequential daynumber during a year. */
int daynumber, /* Sequential daynumber during a year. */
delta_days, /* Whole days since 2000 January 0. */
delta_years; /* Whole years since 2000. */
double cent_J2000, /* Julian centuries since epoch J2000.0 at 0h UT. */
cos_alt, /* Cosine of the altitude of Sun. */
cos_apdec, /* Cosine of the apparent declination of Sun. */
cos_az, /* Cosine of the azimuth of Sun. */
cos_lat, /* Cosine of the site latitude. */
cos_lha, /* Cosine of the local apparent hour angle of Sun. */
days_J2000, /* Days since epoch J2000.0. */
ecliptic_long, /* Solar ecliptic longitude. */
lmst, /* Local mean sidereal time. */
local_ha, /* Local mean hour angle of Sun. */
gmst0h, /* Greenwich mean sidereal time at 0 hours UT. */
integral, /* Integral portion of double precision number. */
mean_anomaly, /* Earth mean anomaly. */
mean_longitude, /* Solar mean longitude. */
mean_obliquity, /* Mean obliquity of the ecliptic. */
pressure = /* Earth mean atmospheric pressure at sea level */
1013.25, /* in millibars. */
sin_apdec, /* Sine of the apparent declination of Sun. */
sin_az, /* Sine of the azimuth of Sun. */
sin_lat, /* Sine of the site latitude. */
tan_alt, /* Tangent of the altitude of Sun. */
temp = /* Earth mean atmospheric temperature at sea level */
15.0, /* in degrees Celsius. */
ut; /* UT hours since midnight. */
/* Check latitude and longitude for proper range before calculating dates.
*/
if (latitude < -90.0 || latitude > 90.0 ||
longitude < -180.0 || longitude > 180.0)
return (-1);
/* If year is not zero then assume date is specified by year, month, day.
* If year is zero then assume date is specified by days_1900.
*/
if (year != 0)
/* Date given by {year, month, day} or {year, 0, daynumber}. */
{
if (year < 1950 || year > 2049)
return (-1);
if (month != 0)
{
if (month < 1 || month > 12 || day < 0.0 || day > 32.0)
return (-1);
daynumber = daynum(year, month, (int)day);
}
else
{
if (day < 0.0 || day > 367.0)
return (-1);
daynumber = (int)day;
}
/* Construct Julian centuries since J2000 at 0 hours UT of date,
* days.fraction since J2000, and UT hours.
*/
delta_years = year - 2000;
/* delta_days is days from 2000/01/00 (1900's are negative). */
delta_days = delta_years * 365 + delta_years / 4 + daynumber;
if (year > 2000)
delta_days += 1;
/* J2000 is 2000/01/01.5 */
days_J2000 = delta_days - 1.5;
cent_J2000 = days_J2000 / 36525.0;
ut = modf(day, &integral);
days_J2000 += ut;
ut *= 24.0;
}
else
/* Date given by days_1900. */
{
/* days_1900 is 18262 for 1950/01/00, and 54788 for 2049/12/32.
* A. A. 1990, K2-K4. */
if (days_1900 < 18262.0 || days_1900 > 54788.0)
return (-1);
/* Construct days.fraction since J2000, UT hours, and
* Julian centuries since J2000 at 0 hours UT of date.
*/
/* days_1900 is 36524 for 2000/01/00. J2000 is 2000/01/01.5 */
days_J2000 = days_1900 - 36525.5;
ut = modf(days_1900, &integral) * 24.0;
cent_J2000 = (integral - 36525.5) / 36525.0;
}
/* Compute solar position parameters.
* A. A. 1990, C24.
*/
mean_anomaly = (357.528 + 0.9856003 * days_J2000);
mean_longitude = (280.460 + 0.9856474 * days_J2000);
/* Put mean_anomaly and mean_longitude in the range 0 -> 2 pi. */
mean_anomaly = modf(mean_anomaly / 360.0, &integral) * TWOPI;
mean_longitude = modf(mean_longitude / 360.0, &integral) * TWOPI;
mean_obliquity = (23.439 - 4.0e-7 * days_J2000) * DEG_RAD;
ecliptic_long = ((1.915 * sin(mean_anomaly)) +
(0.020 * sin(2.0 * mean_anomaly))) * DEG_RAD +
mean_longitude;
*distance = 1.00014 - 0.01671 * cos(mean_anomaly) -
0.00014 * cos(2.0 * mean_anomaly);
/* Tangent of ecliptic_long separated into sine and cosine parts for ap_ra. */
*ap_ra = atan2(cos(mean_obliquity) * sin(ecliptic_long), cos(ecliptic_long));
/* Change range of ap_ra from -pi -> pi to 0 -> 2 pi. */
if (*ap_ra < 0.0)
*ap_ra += TWOPI;
/* Put ap_ra in the range 0 -> 24 hours. */
*ap_ra = modf(*ap_ra / TWOPI, &integral) * 24.0;
*ap_dec = asin(sin(mean_obliquity) * sin(ecliptic_long));
/* Calculate local mean sidereal time.
* A. A. 1990, B6-B7.
*/
/* Horner's method of polynomial exponent expansion used for gmst0h. */
gmst0h = 24110.54841 + cent_J2000 * (8640184.812866 + cent_J2000 *
(0.093104 - cent_J2000 * 6.2e-6));
/* Convert gmst0h from seconds to hours and put in the range 0 -> 24. */
gmst0h = modf(gmst0h / 3600.0 / 24.0, &integral) * 24.0;
if (gmst0h < 0.0)
gmst0h += 24.0;
/* Ratio of lengths of mean solar day to mean sidereal day is 1.00273790934
* in 1990. Change in sidereal day length is < 0.001 second over a century.
* A. A. 1990, B6.
*/
lmst = gmst0h + (ut * 1.00273790934) + longitude / 15.0;
/* Put lmst in the range 0 -> 24 hours. */
lmst = modf(lmst / 24.0, &integral) * 24.0;
if (lmst < 0.0)
lmst += 24.0;
/* Calculate local hour angle, altitude, azimuth, and refraction correction.
* A. A. 1990, B61-B62.
*/
local_ha = lmst - *ap_ra;
/* Put hour angle in the range -12 to 12 hours. */
if (local_ha < -12.0)
local_ha += 24.0;
else if (local_ha > 12.0)
local_ha -= 24.0;
/* Convert latitude and local_ha to radians. */
latitude *= DEG_RAD;
local_ha = local_ha / 24.0 * TWOPI;
cos_apdec = cos(*ap_dec);
sin_apdec = sin(*ap_dec);
cos_lat = cos(latitude);
sin_lat = sin(latitude);
cos_lha = cos(local_ha);
*altitude = asin(sin_apdec * sin_lat + cos_apdec * cos_lha * cos_lat);
cos_alt = cos(*altitude);
/* Avoid tangent overflow at altitudes of +-90 degrees.
* 1.57079615 radians is equal to 89.99999 degrees.
*/
if (fabs(*altitude) < 1.57079615)
tan_alt = tan(*altitude);
else
tan_alt = 6.0e6;
cos_az = (sin_apdec * cos_lat - cos_apdec * cos_lha * sin_lat) / cos_alt;
sin_az = -(cos_apdec * sin(local_ha) / cos_alt);
*azimuth = acos(cos_az);
/* Change range of azimuth from 0 -> pi to 0 -> 2 pi. */
if (atan2(sin_az, cos_az) < 0.0)
*azimuth = TWOPI - *azimuth;
/* Convert ap_dec, altitude, and azimuth to degrees. */
*ap_dec *= RAD_DEG;
*altitude *= RAD_DEG;
*azimuth *= RAD_DEG;
/* Compute refraction correction to be added to altitude to obtain actual
* position.
* Refraction calculated for altitudes of -1 degree or more allows for a
* pressure of 1040 mb and temperature of -22 C. Lower pressure and higher
* temperature combinations yield less than 1 degree refraction.
* NOTE:
* The two equations listed in the A. A. have a crossover altitude of
* 19.225 degrees at standard temperature and pressure. This crossover point
* is used instead of 15 degrees altitude so that refraction is smooth over
* the entire range of altitudes. The maximum residual error introduced by
* this smoothing is 3.6 arc seconds at 15 degrees. Temperature or pressure
* other than standard will shift the crossover altitude and change the error.
*/
/* We want to smooth out the transition for the refraction, rather */
/* than having a discontinuity at zenith = 91 degrees. Thus, we */
/* relax the refraction from it's value at alt == -1 to zero at */
/* alt == -2. trs 3/4/03 */
if (*altitude < -2 || tan_alt == 6.0e6)
*refraction = 0.0;
else if (*altitude < -1)
{
/* 0.241277 * Pres/Temp is refraction at alt == -1 */
/* (alt+2) goes linearly from 1 at alt == -1 to zero */
/* at alt == -2 */
*refraction = 0.241277 * (*altitude + 2) * pressure/(273.0 + temp);
}
else
{
if (*altitude < 19.225)
{
*refraction = (0.1594 + (*altitude) * (0.0196 + 0.00002 * (*altitude))) *
pressure;
*refraction /= (1.0 + (*altitude) * (0.505 + 0.0845 * (*altitude))) *
(273.0 + temp);
}
else
*refraction = 0.00452 * (pressure / (273.0 + temp)) / tan_alt;
}
return 0;
}
/* 'daynum()' returns the sequential daynumber of a calendar date during a
* Gregorian calendar year (for years 1 onward).
* The integer arguments are the four-digit year, the month number, and
* the day of month number.
* (Jan. 1 = 01/01 = 001; Dec. 31 = 12/31 = 365 or 366.)
* A value of -1 is returned if the year is out of bounds.
*/
/* Author: Nels Larson
* Pacific Northwest Lab.
* P.O. Box 999
* Richland, WA 99352
* U.S.A.
*/
int daynum(year, month, day)
int year, month, day;
{
static int begmonth[13] = {0,0,31,59,90,120,151,181,212,243,273,304,334};
int dnum,
leapyr = 0;
/* There is no year 0 in the Gregorian calendar and the leap year cycle
* changes for earlier years. */
if (year < 1)
return (-1);
/* Leap years are divisible by 4, except for centurial years not divisible
* by 400. */
if (((year%4) == 0 && (year%100) != 0) || (year%400) == 0)
leapyr = 1;
dnum = begmonth[month] + day;
if (leapyr && (month > 2))
dnum += 1;
return (dnum);
}
/* "solpos_" is a "fortran jacket" which allows the integer C function
* solarposition() to be called from a fortran program.
*
* Version 1.0 - March 14, 1993.
*
/* Author: Jim Liljegren
* Pacific Northwest Lab.
* P.O. Box 999
* Richland, WA 99352
* U.S.A.
*/
int solpos_(year, month, day, days_1900, latitude, longitude,
ap_ra, ap_dec, altitude, refraction, azimuth, distance)
int *year, /* Four digit year (Gregorian calendar).
* [1950 through 2049; 0 o.k. if using days_1900] */
*month; /* Month number.
* [1 through 12; 0 o.k. if using daynumber for day] */
double *day, /* Calendar day.fraction, or daynumber.fraction.
* [If month is NOT 0:
* 0 through 32; 31st @ 18:10:00 UT = 31.75694
* If month IS 0:
* 0 through 367; 366 @ 18:10:00 UT = 366.75694] */
*days_1900, /* Days since 1900 January 0 @ 00:00:00 UT.
* [18262.0 (1950/01/00) through 54788.0 (2049/12/32);
* 1990/01/01 @ 18:10:00 UT = 32873.75694;
* 0.0 o.k. if using {year, month, day} or
* {year, daynumber}] */
*latitude, /* Observation site geographic latitude.
* [degrees.fraction, North positive] */
*longitude, /* Observation site geographic longitude.
* [degrees.fraction, East positive] */
*ap_ra, /* Apparent solar right ascension.
* [hours; 0.0 <= *ap_ra < 24.0] */
*ap_dec, /* Apparent solar declination.
* [degrees; -90.0 <= *ap_dec <= 90.0] */
*altitude, /* Solar altitude, uncorrected for refraction.
* [degrees; -90.0 <= *altitude <= 90.0] */
*refraction, /* Refraction correction for solar altitude.
* Add this to altitude to compensate for refraction.
* [degrees; 0.0 <= *refraction] */
*azimuth, /* Solar azimuth.
* [degrees; 0.0 <= *azimuth < 360.0, East is 90.0] */
*distance; /* Distance of Sun from Earth (heliocentric-geocentric).
* [astronomical units; 1 a.u. is mean distance] */
{
int retval;
retval= solarposition(*year, *month, *day, *days_1900, *latitude, *longitude,
ap_ra, ap_dec, altitude, refraction, azimuth, distance);
return retval;
}
char *get_solarpos_id()
{
return(rcsid);
}
char *get_solarpos_state()
{
return(rcsstate);
}
src/sunae.c 0 → 100644
View file @ 02287497
#include <math.h>
#include "sunae.h"
#include "la.h"
/*
** Compute solar geometry.
**
** If aep->is_tst is true, then don't do the long
** computation, return GMT;
** Otherwise, do the precise Joe calculation, return TST.
*/
double sunae(ae_pack *aep)
{
double delta, leap, jd, tm, mnlong, mnanom, eclong, oblqec,
num, den, ra, lat, sd, cd, sl, cl, gmst, lmst, refrac,
tst;
int year;
year = aep->year;
if (year < 1950)
year += 1900;
if (aep->is_tst) {
double g, dec, lat, ha, cz, azi, eqtime; /* all RADIANS */
static const double
doy_2_rad = M_2PI/365.0;
lat = M_DTOR*aep->lat;
g = doy_2_rad * (aep->doy - 0.5);
dec = 0.006918 - 0.399912*cos(g) + 0.070257*sin(g) -
0.006758*cos(2.0*g) + 0.000907*sin(2.0*g) -
0.002697*cos(3.0*g) + 0.00148*sin(3.0*g);
aep->dec = M_RTOD*dec;
ha = (0.0833333333*M_PI*(-aep->hour + 12));
aep->ha = M_RTOD * ha;
cz = sin(dec)*sin(lat) + cos(dec)*cos(lat)*cos(ha);
aep->el = M_RTOD * (M_PI_2 - acos(cz));
azi = (cz*sin(lat) - sin(dec)) / (sin(acos(cz))*cos(lat));
aep->az = M_RTOD*azi;
eqtime = 0.000075 + 0.001868*cos(g) - 0.032077*sin(g)
- 0.014615*cos(2.0*g) - 0.040848*sin(2.0*g);
return aep->hour - M_RTOH*eqtime;
}
/*
** the time is gmt: do the hard calculation
*/
delta = year - 1949.0;
leap = (int) (0.25 * delta);
jd = 32916.5 + delta*365.0 + leap + aep->doy + aep->hour/24.0;
/*
** "the last year of a century is not a leap year therefore"
*/
if ((year % 100 == 0) ) /* && (year % 4 != 0)) */
jd -= 1.0;
tm = jd - 51545.0;
mnlong = fmod(280.460 + 0.9856474*tm, 360.0);
if (mnlong < 0.0)
mnlong += 360.0;
mnanom = fmod(357.528 + 0.9856003*tm, 360.0);
mnanom *= M_DTOR;
eclong = fmod(mnlong + 1.915*sin(mnanom) + 0.020*sin(2.0*mnanom), 360.0);
eclong *= M_DTOR;
oblqec = 23.439 - 0.0000004*tm;
oblqec *= M_DTOR;
num = cos(oblqec)*sin(eclong);
den = cos(eclong);
ra = atan(num/den);
/* force ra between 0 and 2*pi */
if (den < 0.0)
ra += M_PI;
else if (num < 0.0)
ra += M_2PI;
/* dec in radians */
aep->dec = asin(sin(oblqec)*sin(eclong));
/* calculate Greenwich mean sidereal time in hours */
gmst = fmod(6.697375 + 0.0657098242*tm + aep->hour, 24.0);
/* calculate local mean sidereal time in radians */
lmst = fmod(gmst + aep->lon*M_DTOH, 24.0);
lmst = lmst*M_HTOR;
/* calculate hour angle in radians between -pi and pi */
aep->ha = lmst - ra;
if (aep->ha < -M_PI)
aep->ha += M_2PI;
if (aep->ha > M_PI)
aep->ha -= M_2PI;
/* change latitude to radians */
lat = aep->lat*M_DTOR;
sd = sin(aep->dec);
cd = cos(aep->dec);
sl = sin(lat);
cl = cos(lat);
/* calculate azimuth and elevation */
aep->el = asin(sd*sl + cd*cl*cos(aep->ha));
aep->az = asin(-cd*sin(aep->ha)/cos(aep->el));
/* this puts azimuth between 0 and 2*pi radians */
if (sd - sin(aep->el)*sl >= 0.0) {
if(sin(aep->az) < 0.0)
aep->az += M_2PI;
} else
aep->az = M_PI - aep->az;
/* calculate refraction correction for US stan. atmosphere */
/* need to have el in degs before calculating correction */
aep->el *= M_RTOD;
if (aep->el >= 19.225)
refrac = 0.00452*3.51823/tan(M_DTOR*aep->el);
else if (aep->el > -0.766 && aep->el < 19.225)
refrac = 3.51823*(0.1594 + 0.0196*aep->el + 0.00002*aep->el*aep->el) /
(1.0 + 0.505*aep->el + 0.0845*aep->el*aep->el);
else
refrac = 0.0;
/*
if (aep->el > -0.762)
refrac = 3.51561*(0.1594 + 0.0196*aep->el + 0.00002*aep->el*aep->el) /
(1.0 + 0.505*aep->el + 0.0845*aep->el*aep->el);
else
refrac = 0.762;
*/
/* note that 3.51561=1013.2 mb/288.2 C */
aep->el += refrac;
/* calculate distance to sun in A.U. & diameter in degs */
aep->soldst = 1.00014 - 0.01671*cos(mnanom) - 0.00014*cos(mnanom + mnanom);
/* convert az and lat to degs before returning */
aep->az *= M_RTOD;
aep->ha *= M_RTOD;
aep->dec *= M_RTOD;
tst = 12.0 + aep->ha*M_DTOH;
while (tst < 0.0)
tst += 24.0;
while (tst > 24.0)
tst -= 24.0;
return tst;
}
/*
** this function calculates air mass using kasten's
** approximation to bemporad's tables
*/
double armass(double el)
{
double zr;
zr = M_DTOR * (90.0 - el);
return 1.0/(cos(zr) + 0.50572*pow(6.07995 + el, -1.6364));
}
src/sunae.h 0 → 100644
View file @ 02287497
#ifndef AE_H
#define AE_H
typedef struct {
int year, doy, is_tst;
double hour, lat, lon,
az, el, ha, dec, soldst;
} ae_pack;
double sunae(ae_pack *),
armass(double);
#endif /* AE_H */
src/time_of.c 0 → 100644
View file @ 02287497
#include <stdio.h>
#include "la.h"
#include "util.h"
/*
** return the time_t corresponding to the pointed-to obs_time_t
*/
time_t time_of(obs_time_t *r)
{
struct tm stm, *dur;
time_t t,s,gmktime();
long daysec;
int mo, dom;
doy2md(r->year, r->doy, &mo, &dom);
stm.tm_isdst = 0;
stm.tm_year = r->year;
stm.tm_hour = (int) r->gmt;
daysec = (long) (3600*r->gmt);
stm.tm_min = (daysec % 3600)/60;
stm.tm_sec = daysec % 60;
stm.tm_mday = dom;
stm.tm_mon = mo - 1;
/* t = mktime(&stm); */ /* mktime() interprets stm in local time */
t = gmktime(&stm);
/* dur = gmtime(&t); */ /* debugger test */
#ifdef CHECK_DATE
if (stm.tm_yday + 1 != r->doy) { /* just a sanity check */
fprintf(stderr, "time_of: doy problem; %d - %d\n", stm.tm_yday,
r->doy);
exit(1);
}
#endif
return t;
}
static char reg_days[] = { 31, 28, 31, 30, 31, 30,
31, 31, 30, 31, 30, 31 };
static char leap_days[] = { 31, 29, 31, 30, 31, 30,
31, 31, 30, 31, 30, 31 };
#define HRS_PER_YEAR 24 * 365L
#define DAYS_1900_1970 25568L /* days from 1/1/00 to 1/1/1970 */
/*
* gmktime() emulates the mktime() function but whereas mktime is intended to be
* a sort of inverse of localtime, gmktime is an inverse of gmtime. Called the same
* way, though.
*/
time_t gmktime(struct tm *stm)
{
char *days_per_mon;
time_t jdays, yr, secs=0;
int i;
/* vulnerable to the dreaded 2071 bug! */
if(stm->tm_year < 70)
yr = stm->tm_year + 2000;
else
yr = stm->tm_year + 1900;
/* this def of yeap year will be good for a long time */
if(yr % 4 == 0)
days_per_mon = leap_days;
else
days_per_mon = reg_days;
jdays = stm->tm_mday; /* start by taking the number of days this month */
for(i=0;i<stm->tm_mon;i++)
jdays += days_per_mon[i]; /* add in all the previous months */
jdays += ((long) yr - 1900) * 365; /* days per year * years since 1900 */
jdays += ((long) yr - 1900 - 1)/4; /* leap days */
jdays -= (time_t) DAYS_1900_1970; /* reference to "Unix start time" */
secs = jdays * 3600 * 24; /* crank out the seconds */
secs += stm->tm_sec + stm->tm_min*60 + stm->tm_hour*3600; /* add in today's time */
return(secs);
}
src/util.c 0 → 100644
View file @ 02287497
int md2doy(int yr, int mon, int dom)
{
static const int dim_n[12] = {31, 28, 31, 30, 31, 30,
31, 31, 30, 31, 30, 31},
dim_l[12] = {31, 29, 31, 30, 31, 30,
31, 31, 30, 31, 30, 31};
int i, *dim, doy;
if (mon < 1 || mon > 12 || dom < 1 || dom > 31)
return 0;
if ((yr % 4 == 0 && yr % 100 != 0) || yr %400 == 0)
dim = (int *) dim_l;
else
dim = (int *) dim_n;
for (doy=0, i=0; mon - 1 > i; i++)
doy += dim[i];
if (dom > dim[i])
return 0;
doy += dom;
return doy;
}
void doy2md(int yr, int doy, int *mo, int *dom)
{
static const int ldm_n[14] = {0, 0, 31, 59, 90, 120, 151, 181,
212, 243, 273, 304, 334, 365},
ldm_l[14] = {0, 0, 31, 60, 91, 121, 152, 182,
213, 244, 274, 305, 335, 366};
int *ldm;
if (doy < 1 || doy > 366) {
*mo = 0;
*dom = 0;
return;
}
yr += 1900; /* yr is 100 in 2000, as in struct tm */
if ((yr % 4 == 0 && yr % 100 != 0) || yr %400 == 0)
ldm = (int *) ldm_l;
else {
if (doy == 366) {
*mo = 0;
*dom = 0;
return;
}
ldm = (int *) ldm_n;
}
for (*mo=1; doy > ldm[*mo + 1]; (*mo)++) {
}
*dom = doy - ldm[*mo];
}
/*
** split - divide a string into fields, like awk split()
** return number of fields, including overflow
** sep: "" white, "c" single char, "ab" [ab]+
*/
int split(char *string, char *fields[], int nfields, char *sep)
{
char *p = string, c, sepc2, *sepp,
sepc = sep[0],
**fp = fields;
int fn, trimtrail;
/* white space */
if (sepc == '\0') {
while ((c = *p++) == ' ' || c == '\t')
continue;
p--;
trimtrail = 1;
sep = " \t"; /* note, code below knows this is 2 long */
sepc = ' ';
} else
trimtrail = 0;
sepc2 = sep[1]; /* now we can safely pick this up */
/* catch empties */
if (*p == '\0')
return(0);
/* single separator */
if (sepc2 == '\0') {
fn = nfields;
for (;;) {
*fp++ = p;
fn--;
if (fn == 0)
break;
while ((c = *p++) != sepc)
if (c == '\0')
return(nfields - fn);
*(p-1) = '\0';
}
/* we have overflowed the fields vector -- just count them */
fn = nfields;
for (;;) {
while ((c = *p++) != sepc)
if (c == '\0')
return(fn);
fn++;
}
/* not reached */
}
/* two separators */
if (sep[2] == '\0') {
fn = nfields;
for (;;) {
*fp++ = p;
fn--;
while ((c = *p++) != sepc && c != sepc2)
if (c == '\0') {
if (trimtrail && **(fp-1) == '\0')
fn++;
return(nfields - fn);
}
if (fn == 0)
break;
*(p-1) = '\0';
while ((c = *p++) == sepc || c == sepc2)
continue;
p--;
}
/* we have overflowed the fields vector -- just count them */
fn = nfields;
while (c != '\0') {
while ((c = *p++) == sepc || c == sepc2)
continue;
p--;
fn++;
while ((c = *p++) != '\0' && c != sepc && c != sepc2)
continue;
}
/* might have to trim trailing white space */
if (trimtrail) {
p--;
while ((c = *--p) == sepc || c == sepc2)
continue;
p++;
if (*p != '\0') {
if (fn == nfields+1)
*p = '\0';
fn--;
}
}
return(fn);
}
/* n separators */
fn = 0;
for (;;) {
if (fn < nfields)
*fp++ = p;
fn++;
for (;;) {
c = *p++;
if (c == '\0')
return(fn);
sepp = sep;
while ((sepc = *sepp++) != '\0' && sepc != c)
continue;
if (sepc != '\0') /* it was a separator */
break;
}
if (fn < nfields)
*(p-1) = '\0';
for (;;) {
c = *p++;
sepp = sep;
while ((sepc = *sepp++) != '\0' && sepc != c)
continue;
if (sepc == '\0') /* it wasn't a separator */
break;
}
p--;
}
/* not reached */
}
src/util.h 0 → 100644
View file @ 02287497
#ifndef UTIL_H
#define UTIL_H
int split(char *string, char *fields[], int nfields, char *sep);
void doy2md(int y, int doy, int *mo, int *dom);
#endif
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment