|
|
/* $Id: geometry.c 4554 2014-12-17 15:47:14Z brideout $ */
|
|
|
|
|
|
#include <math.h>
|
|
|
#include <stdlib.h>
|
|
|
#include <string.h>
|
|
|
#include <sys/types.h>
|
|
|
#include <ctype.h>
|
|
|
#include <sys/stat.h>
|
|
|
#include <fcntl.h>
|
|
|
#include <stdio.h>
|
|
|
#include <cedarIO.h>
|
|
|
#include <madrec.h>
|
|
|
#include <cedar.h>
|
|
|
#include <date.h>
|
|
|
|
|
|
/* include code written by National Renewable Energy Lab */
|
|
|
#include <solpos00.h>
|
|
|
|
|
|
/***********************************************************************
|
|
|
*
|
|
|
* sprod calculates the scalar product of two vectors a and b,
|
|
|
* sprod = a .dot. b.
|
|
|
*/
|
|
|
double
|
|
|
sprod(double *a, double *b)
|
|
|
{
|
|
|
double sp;
|
|
|
|
|
|
sp = a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
|
|
|
return(sp);
|
|
|
}
|
|
|
|
|
|
|
|
|
/***********************************************************************
|
|
|
*
|
|
|
* vadd calculates the sum of two vectors a and b, c = a + b.
|
|
|
*/
|
|
|
int
|
|
|
vadd(double *a, double *b, double *c)
|
|
|
{
|
|
|
c[0] = a[0] + b[0];
|
|
|
c[1] = a[1] + b[1];
|
|
|
c[2] = a[2] + b[2];
|
|
|
return (0);
|
|
|
}
|
|
|
|
|
|
|
|
|
/***********************************************************************
|
|
|
*
|
|
|
* vsub calculates the difference of two vectors a and b, c = a - b.
|
|
|
*/
|
|
|
int
|
|
|
vsub(double *a, double *b, double *c)
|
|
|
{
|
|
|
c[0] = a[0] - b[0];
|
|
|
c[1] = a[1] - b[1];
|
|
|
c[2] = a[2] - b[2];
|
|
|
return(0);
|
|
|
}
|
|
|
|
|
|
|
|
|
/***********************************************************************
|
|
|
*
|
|
|
* csconv converts between cartesian coordinates x,y,z and spherical
|
|
|
* coordinates r,theta,phi. if imode=1, (x,y,z) -> (r,theta,phi).
|
|
|
* if imode=2, (r,theta,phi) -> (x,y,z). theta and phi are in
|
|
|
* degrees.
|
|
|
*/
|
|
|
int
|
|
|
csconv(double *xp, double *yp, double *zp,
|
|
|
double *rp, double *thetap, double *phip,
|
|
|
int imode)
|
|
|
{
|
|
|
double dtr=0.0174532925199, rlmin=1.e-20,
|
|
|
x, y, z, r, theta, phi,
|
|
|
rho2, ct, t, st, cp, p, sp;
|
|
|
|
|
|
if (imode == 1) {
|
|
|
x = *xp;
|
|
|
y = *yp;
|
|
|
z = *zp;
|
|
|
rho2 = x*x + y*y;
|
|
|
r = sqrt(rho2 + z*z);
|
|
|
if (z < rlmin && z > 0)
|
|
|
z = rlmin;
|
|
|
else if (z > -rlmin && z < 0)
|
|
|
z = -rlmin;
|
|
|
else
|
|
|
z = z;
|
|
|
theta = atan2(sqrt(rho2), z)/dtr;
|
|
|
if (x < rlmin && x > 0)
|
|
|
x = rlmin;
|
|
|
else if (x > -rlmin && x < 0)
|
|
|
x = -rlmin;
|
|
|
else
|
|
|
x = x;
|
|
|
phi = atan2(y, x)/dtr;
|
|
|
*rp = r;
|
|
|
*thetap = theta;
|
|
|
*phip = phi;
|
|
|
return(0);
|
|
|
}
|
|
|
else if (imode == 2) {
|
|
|
r = *rp;
|
|
|
theta = *thetap;
|
|
|
phi = *phip;
|
|
|
t = dtr*theta;
|
|
|
p = dtr*phi;
|
|
|
ct = cos(t);
|
|
|
st = sin(t);
|
|
|
cp = cos(p);
|
|
|
sp = sin(p);
|
|
|
x = r*st*cp;
|
|
|
y = r*st*sp;
|
|
|
z = r*ct;
|
|
|
*xp = x;
|
|
|
*yp = y;
|
|
|
*zp = z;
|
|
|
return(0);
|
|
|
}
|
|
|
else {
|
|
|
return(1);
|
|
|
}
|
|
|
}
|
|
|
|
|
|
/***********************************************************************
|
|
|
*
|
|
|
* vctcnv converts between the cartesian and spherical coordinate
|
|
|
* representations of a vector field f. (fx,fy,fz) are the
|
|
|
* components of the field at (x,y,z). (fr,ft,fp) are the
|
|
|
* components of the field at (r,theta,phi) in the directions of
|
|
|
* increasing r, increasing theta and increasing phi. if imode=1,
|
|
|
* (fx,fy,fz,x,y,z) -> (fr,ft,fp,r,theta,phi). if imode=2,
|
|
|
* (fr,ft,fp,r,theta,phi) -> (fx,fy,fz,x,y,z). theta and phi are
|
|
|
* in degrees.
|
|
|
*/
|
|
|
int
|
|
|
vctcnv(double *fxp, double *fyp, double *fzp,
|
|
|
double *xp, double *yp, double *zp,
|
|
|
double *frp, double *ftp, double *fpp,
|
|
|
double *rp, double *thetap, double *phip,
|
|
|
int imode)
|
|
|
{
|
|
|
double dtr=0.0174532925199, rlmin=1.e-20,
|
|
|
fx, fy, fz, x, y, z, fr, ft, fp, r, theta, phi,
|
|
|
rho2, ct, st, cp, sp, t, p;
|
|
|
|
|
|
if (imode == 1) {
|
|
|
fx = *fxp;
|
|
|
fy = *fyp;
|
|
|
fz = *fzp;
|
|
|
x = *xp;
|
|
|
y = *yp;
|
|
|
z = *zp;
|
|
|
rho2 = x*x + y*y;
|
|
|
r = sqrt(rho2 + z*z);
|
|
|
if (z < rlmin && z > 0)
|
|
|
z = rlmin;
|
|
|
else if (z > -rlmin && z < 0)
|
|
|
z = -rlmin;
|
|
|
else
|
|
|
z = z;
|
|
|
theta = atan2(sqrt(rho2), z);
|
|
|
if (x < rlmin && x > 0)
|
|
|
x = rlmin;
|
|
|
else if (x > -rlmin && x < 0)
|
|
|
x = -rlmin;
|
|
|
else
|
|
|
x = x;
|
|
|
phi = atan2(y, x);
|
|
|
ct = cos(theta);
|
|
|
st = sin(theta);
|
|
|
cp = cos(phi);
|
|
|
sp = sin(phi);
|
|
|
theta = theta/dtr;
|
|
|
phi = phi/dtr;
|
|
|
fr = st*cp*fx + st*sp*fy + ct*fz;
|
|
|
ft = ct*cp*fx + ct*sp*fy - st*fz;
|
|
|
fp = -sp*fx + cp*fy;
|
|
|
*frp= fr;
|
|
|
*ftp = ft;
|
|
|
*fpp = fp;
|
|
|
*rp= r;
|
|
|
*thetap = theta;
|
|
|
*phip = phi;
|
|
|
return(0);
|
|
|
}
|
|
|
else if (imode == 2) {
|
|
|
fr = *frp;
|
|
|
ft = *ftp;
|
|
|
fp = *fpp;
|
|
|
r = *rp;
|
|
|
theta = *thetap;
|
|
|
phi = *phip;
|
|
|
t = dtr*theta;
|
|
|
p = dtr*phi;
|
|
|
ct = cos(t);
|
|
|
st = sin(t);
|
|
|
cp = cos(p);
|
|
|
sp = sin(p);
|
|
|
x = r*st*cp;
|
|
|
y = r*st*sp;
|
|
|
z = r*ct;
|
|
|
fx = st*cp*fr + ct*cp*ft - sp*fp;
|
|
|
fy = st*sp*fr + ct*sp*ft + cp*fp;
|
|
|
fz = ct*fr - st*ft;
|
|
|
*fxp = fx;
|
|
|
*fyp = fy;
|
|
|
*fzp = fz;
|
|
|
*xp = x;
|
|
|
*yp = y;
|
|
|
*zp = z;
|
|
|
return(0);
|
|
|
}
|
|
|
else {
|
|
|
return(1);
|
|
|
}
|
|
|
}
|
|
|
|
|
|
/***********************************************************************
|
|
|
*
|
|
|
* point calculates the position of a point defined by the radar
|
|
|
* line-of sight vector to that point.
|
|
|
*
|
|
|
* input parameters
|
|
|
* sr - distance of station from center of earth (km)
|
|
|
* slat - geocentric latitude of station (deg)
|
|
|
* slon - longitude of station (deg)
|
|
|
* az - radar azimuth (deg)
|
|
|
* el - radar elevation (deg)
|
|
|
* range - radar range (km)
|
|
|
*
|
|
|
* output parameters
|
|
|
* pr - distance from center of earth of observation point (km)
|
|
|
* glat - observation point geocentric latitude (deg)
|
|
|
* glon - observation point longitude (deg)
|
|
|
*/
|
|
|
int
|
|
|
point(double *srp, double *slatp, double *slonp,
|
|
|
double *azp, double *elp, double *rangep,
|
|
|
double *prp, double *glatp, double *glonp)
|
|
|
{
|
|
|
double sr, slat, slon, az, el, range, pr, glat, glon,
|
|
|
s[3], r[3], p[3], rt, rp, rr, t, el1, az1, slat1;
|
|
|
|
|
|
sr = *srp;
|
|
|
slat = *slatp;
|
|
|
slon = *slonp;
|
|
|
az = *azp;
|
|
|
el = *elp;
|
|
|
range = *rangep;
|
|
|
|
|
|
/* calculate "line-of-sight" station centered cartesian coords */
|
|
|
el1 = 90.0 - el;
|
|
|
az1 = 180.0 - az;
|
|
|
(void) csconv(&rt, &rp, &rr, &range, &el1, &az1, 2);
|
|
|
|
|
|
/* calculate "line-of-sight" earth centered cartesian coords
|
|
|
and "station" earth centered cartesian coords */
|
|
|
slat1 = 90.0 - slat;
|
|
|
(void) vctcnv(&r[0], &r[1], &r[2], &s[0], &s[1], &s[2],
|
|
|
&rr, &rt, &rp, &sr, &slat1, &slon, 2);
|
|
|
|
|
|
/* calculate "observation-point" earth centered cartesian coords */
|
|
|
(void) vadd(s, r, p);
|
|
|
|
|
|
/* calculate "observation-point" earth centered spherical coords */
|
|
|
(void) csconv(&p[0], &p[1], &p[2], &pr, &t, &glon, 1);
|
|
|
glat = 90. - t;
|
|
|
|
|
|
*prp = pr;
|
|
|
*glatp = glat;
|
|
|
*glonp = glon;
|
|
|
return(0);
|
|
|
}
|
|
|
|
|
|
|
|
|
/***********************************************************************
|
|
|
*
|
|
|
* look calculates the azimuth, elevation and range from a radar
|
|
|
* of a specified point.
|
|
|
*
|
|
|
* input parameters
|
|
|
* sr - distance of station from center of earth (km)
|
|
|
* slat - geocentric latitude of station (deg)
|
|
|
* slon - longitude of station (deg)
|
|
|
* pr - distance from center of earth of observation point (km)
|
|
|
* glat - observation point geocentric latitude (deg)
|
|
|
* glon - observation point longitude (deg)
|
|
|
*
|
|
|
* output parameters
|
|
|
* az - radar azimuth (deg)
|
|
|
* el - radar elevation (deg)
|
|
|
* range - radar range (km)
|
|
|
*/
|
|
|
int
|
|
|
look(double *srp, double *slatp, double *slonp,
|
|
|
double *prp, double *glatp, double *glonp,
|
|
|
double *azp, double *elp, double *rangep)
|
|
|
{
|
|
|
double sr, slat, slon, pr, glat, glon, az, el, range,
|
|
|
s[3], r[3], p[3], rr, rt, rp, sr1, st2, sp1, glat1, slat1;
|
|
|
|
|
|
sr = *srp;
|
|
|
slat = *slatp;
|
|
|
slon = *slonp;
|
|
|
pr = *prp;
|
|
|
glat = *glatp;
|
|
|
glon = *glonp;
|
|
|
|
|
|
/*calculate "observation-point" earth centered cartesian coords */
|
|
|
glat1 = 90.0 - glat;
|
|
|
(void) csconv(&p[0], &p[1], &p[2], &pr, &glat1, &glon, 2);
|
|
|
|
|
|
/*calculate "station" earth centered cartesian coordinates*/
|
|
|
slat1 = 90.0 - slat;
|
|
|
(void) csconv(&s[0], &s[1], &s[2], &sr, &slat1, &slon, 2);
|
|
|
|
|
|
/*calculate "line-of-sight" earth centered cartesian coords*/
|
|
|
(void) vsub(p, s, r);
|
|
|
|
|
|
/*calculate "line-of-sight" station centered cartesian coords*/
|
|
|
(void) vctcnv(&r[0], &r[1], &r[2], &s[0], &s[1], &s[2],
|
|
|
&rr, &rt, &rp, &sr1, &st2, &sp1, 1);
|
|
|
|
|
|
/*calculate "line-of-sight" station centered spherical coords*/
|
|
|
(void) csconv(&rt, &rp, &rr, &range, &el, &az, 1);
|
|
|
el = 90. - el;
|
|
|
az = 180. - az;
|
|
|
|
|
|
*azp = az;
|
|
|
*elp = el;
|
|
|
*rangep = range;
|
|
|
return(0);
|
|
|
}
|
|
|
|
|
|
|
|
|
/***********************************************************************
|
|
|
*
|
|
|
* convrt converts between geodetic and geocentric coordinates. the
|
|
|
* reference geoid is that adopted by the iau in 1964. a=6378.16,
|
|
|
* b=6356.7746, f=1/298.25. the equations for conversion from
|
|
|
* geocentric to geodetic are from astron. j., vol 66, 1961, p. 15.
|
|
|
* i=1 geodetic to geocentric
|
|
|
* i=2 geocentric to geodetic
|
|
|
* gdlat geodetic latitude (degrees)
|
|
|
* gdalt altitude above geoid (km)
|
|
|
* gclat geocentric latitude (degrees)
|
|
|
* rkm geocentric radial distance (km)
|
|
|
*/
|
|
|
int
|
|
|
convrt(int i, double *gdlatp, double *gdaltp,
|
|
|
double *gclatp, double *rkmp)
|
|
|
|
|
|
{
|
|
|
double gdlat, gdalt, gclat, rkm,
|
|
|
sinlat,gdl,coslat,cl2,sb2,sinbet,cosbet,
|
|
|
x,y,a2,rer,ccl,gcl,scl,s4cl,s2cl,c2cl,sl2,rgeoid,a4,a6,
|
|
|
a8,c4cl,s8cl,s6cl,dltcl;
|
|
|
double a=6378.16, ab2=1.0067397, ep2=0.0067397, dtr=0.0174532925199;
|
|
|
|
|
|
/* geodetic to geocentric */
|
|
|
if (i == 1) {
|
|
|
gdlat = *gdlatp;
|
|
|
gdalt = *gdaltp;
|
|
|
gdl = dtr*gdlat;
|
|
|
sinlat = sin(gdl);
|
|
|
coslat = cos(gdl);
|
|
|
sl2 = sinlat*sinlat;
|
|
|
cl2 = ab2*coslat;
|
|
|
cl2 = cl2*cl2;
|
|
|
sinbet = sinlat/sqrt(sl2+cl2);
|
|
|
sb2 = sinbet*sinbet;
|
|
|
if (sb2 > 1.0) sb2=1.0;
|
|
|
cosbet = sqrt(1. - sb2);
|
|
|
rgeoid = a/sqrt(1. + ep2*sb2);
|
|
|
x = rgeoid*cosbet + gdalt*coslat;
|
|
|
y = rgeoid*sinbet + gdalt*sinlat;
|
|
|
rkm = sqrt(x*x + y*y);
|
|
|
gclat = atan2(y,x)/dtr;
|
|
|
*gclatp = gclat;
|
|
|
*rkmp = rkm;
|
|
|
return(0);
|
|
|
}
|
|
|
/* geocentic to geodetic */
|
|
|
else if (i == 2) {
|
|
|
gclat = *gclatp;
|
|
|
rkm = *rkmp;
|
|
|
rer = rkm/a;
|
|
|
a2 = ((-1.4127348e-8/rer + .94339131e-8)/rer + .33523288e-2)/rer;
|
|
|
a4 = (((-1.2545063e-10/rer + .11760996e-9)/rer +
|
|
|
.11238084e-4)/rer - .2814244e-5)/rer;
|
|
|
a6 = ((54.939685e-9/rer - 28.301730e-9)/rer + 3.5435979e-9)/rer;
|
|
|
a8 = (((320./rer - 252.)/rer + 64.)/rer - 5.)/rer*.98008304e-12;
|
|
|
gcl = dtr*gclat;
|
|
|
ccl = cos(gcl);
|
|
|
scl = sin(gcl);
|
|
|
s2cl = 2.*scl*ccl;
|
|
|
c2cl = 2.*ccl*ccl - 1.0;
|
|
|
s4cl = 2.*s2cl*c2cl;
|
|
|
c4cl = 2.*c2cl*c2cl - 1.0;
|
|
|
s8cl = 2.*s4cl*c4cl;
|
|
|
s6cl = s2cl*c4cl + c2cl*s4cl;
|
|
|
dltcl = s2cl*a2 + s4cl*a4 + s6cl*a6 + s8cl*a8;
|
|
|
gdlat = gclat + dltcl/dtr;
|
|
|
gdalt = rkm - a/sqrt(1.+ep2*scl*scl);
|
|
|
*gdlatp = gdlat;
|
|
|
*gdaltp = gdalt;
|
|
|
return(0);
|
|
|
}
|
|
|
else {
|
|
|
return(1);
|
|
|
}
|
|
|
}
|
|
|
|
|
|
|
|
|
/***********************************************************************
|
|
|
*
|
|
|
* rpcart computes the components (rfx,rfy,rfz) relative to an earth
|
|
|
* centered cartesian coordinate system of the radar line of sight
|
|
|
* vector from a radar with coordinates sr (distance from center
|
|
|
* of earth), slat (geocentric latitude) and slon (longitude). the
|
|
|
* observation point is specified by az (azimuth), el (elevation) and
|
|
|
* range (range). the cartesian coordinates of the observation
|
|
|
* point are returned in (pfx,pfy,pfz).
|
|
|
* input - sr,slat,slon,az,el,range
|
|
|
* output - rfx,rfy,rfz,pfx,pfy,pfz
|
|
|
*/
|
|
|
int
|
|
|
rpcart (double *srp, double *slatp, double *slonp,
|
|
|
double *azp, double *elp, double *rangep,
|
|
|
double *rfxp, double *rfyp, double *rfzp,
|
|
|
double *pfxp, double *pfyp, double *pfzp)
|
|
|
{
|
|
|
double dtr=0.0174532925199, rr, rtheta, rphi, a, e,
|
|
|
ca, sa, ce, se, rx, ry, rz, rfr, rft, rfp;
|
|
|
|
|
|
rr = *srp;
|
|
|
rtheta = 90. - *slatp;
|
|
|
rphi = *slonp;
|
|
|
a = dtr*(180. - *azp);
|
|
|
e = dtr*(90. - *elp);
|
|
|
ca = cos(a);
|
|
|
sa = sin(a);
|
|
|
ce = cos(e);
|
|
|
se = sin(e);
|
|
|
rfr = *rangep*ce;
|
|
|
rft = *rangep*se*ca;
|
|
|
rfp = *rangep*se*sa;
|
|
|
(void) vctcnv(rfxp, rfyp, rfzp, &rx, &ry, &rz, &rfr, &rft, &rfp, &rr,
|
|
|
&rtheta, &rphi, 2);
|
|
|
*pfxp = rx + *rfxp;
|
|
|
*pfyp = ry + *rfyp;
|
|
|
*pfzp = rz + *rfzp;
|
|
|
return(0);
|
|
|
}
|
|
|
|
|
|
|
|
|
/***********************************************************************
|
|
|
*
|
|
|
* gdv converts a vector field f at geodetic latitude gdlat and
|
|
|
* geocentric latitude gclat from a geocentric based representation
|
|
|
* to a geodetic based representation. the geocentric components
|
|
|
* are fr (radial outward), ft (increasing geocentric colatitude,
|
|
|
* e.g. southward) and fp (increasing east longitude). the
|
|
|
* geodetic components are fx (northward, parallel to surface of
|
|
|
* earth), fy (eastward, parallel to surface of earth) and fz
|
|
|
* (downward, perpendicular to surface of earth). fr,ft,fp thus
|
|
|
* correspond to spherical coordinates r,theta,phi, with their
|
|
|
* origin at the center of the earth. x,y,z are the coordinates
|
|
|
* customarily used to describe the three components of the
|
|
|
* geomagnetic field. fp and fy are the same.
|
|
|
*/
|
|
|
int
|
|
|
gdv(double *gdlatp, double *gclatp,
|
|
|
double *frp, double *ftp, double *fpp,
|
|
|
double *fxp, double *fyp, double *fzp)
|
|
|
{
|
|
|
double dtr=0.0174532925199, gdl, sinlat, coslat, t,
|
|
|
ct, st, sind, cosd;
|
|
|
|
|
|
gdl = dtr*(*gdlatp);
|
|
|
sinlat = sin(gdl);
|
|
|
coslat = cos(gdl);
|
|
|
t = dtr*(90. - *gclatp);
|
|
|
ct = cos(t);
|
|
|
st = sin(t);
|
|
|
sind = st*sinlat - ct*coslat;
|
|
|
cosd = ct*sinlat + st*coslat;
|
|
|
*fxp = -*ftp*cosd - *frp*sind;
|
|
|
*fyp = *fpp;
|
|
|
*fzp = *ftp*sind - *frp*cosd;
|
|
|
return (0);
|
|
|
}
|
|
|
|
|
|
|
|
|
/***********************************************************************
|
|
|
*
|
|
|
* los2geodetic calculates the position of a point defined by an instrument
|
|
|
* line-of sight vector to that point. This is a convenience routine in
|
|
|
* which the instrument location is specified by its CEDAR instrument code
|
|
|
* and which returns the geodetic coordinates of the point.
|
|
|
*
|
|
|
*
|
|
|
* input parameters
|
|
|
* kinst - instrument code in metadata
|
|
|
* az - radar azimuth (deg)
|
|
|
* el - radar elevation (deg)
|
|
|
* range - radar range (km)
|
|
|
*
|
|
|
* output parameters
|
|
|
* gdlat - observation point geodetic latitude (deg)
|
|
|
* glon - observation point longitude (deg)
|
|
|
* gdalt - altitude above geoid (km)
|
|
|
*/
|
|
|
|
|
|
int los2geodetic(int kinst, double az, double el, double range,
|
|
|
double *gdlatp, double *glonp, double *gdaltp)
|
|
|
{
|
|
|
double slat, slon, salt, sgclat, sr, pr, pgclat, plat, plon, palt;
|
|
|
|
|
|
cedarGetStationPos(kinst, &slat, &slon, &salt);
|
|
|
|
|
|
/* check for missing info (salt never missing) */
|
|
|
if (isnan(slat) || isnan(salt)) {
|
|
|
*gdlatp = missing;
|
|
|
*glonp = missing;
|
|
|
*gdaltp = missing;
|
|
|
return(-1);
|
|
|
}
|
|
|
|
|
|
/* Convert station coordinates from geodetic to geocentric */
|
|
|
(void) convrt(1, &slat, &salt, &sgclat, &sr);
|
|
|
|
|
|
/* Compute geocentric coordinates of specified point */
|
|
|
(void) point(&sr, &sgclat, &slon, &az, &el, &range,
|
|
|
&pr, &pgclat, &plon);
|
|
|
|
|
|
/* Convert observation point coordinates from geocentric to geodetic */
|
|
|
(void) convrt(2, &plat, &palt, &pgclat, &pr);
|
|
|
|
|
|
*gdlatp = plat;
|
|
|
*glonp = plon;
|
|
|
*gdaltp = palt;
|
|
|
|
|
|
return(0);
|
|
|
}
|
|
|
|
|
|
|
|
|
/***********************************************************************
|
|
|
*
|
|
|
* solarzen_az calculates the solar zenith and az angles for a given time, gdlat,
|
|
|
* and glon.
|
|
|
*
|
|
|
*
|
|
|
* input parameters
|
|
|
* double ut - Universal time in seconds since 1950
|
|
|
* double gdlat - geodetic latitude in degrees
|
|
|
* double glon - geodetic longitude in degrees
|
|
|
*
|
|
|
* output parameters
|
|
|
* szen - solar zenith angle (deg, 0=directly overhead)
|
|
|
* saz - solar azimuth angle (deg, N=0, E=90)
|
|
|
*
|
|
|
*
|
|
|
* Solar zenith angle is calculated at 0 alt, although this changes
|
|
|
* very little with altitute. No atmospheric correction is applied.
|
|
|
* This method uses solpos.c, written by National Renewable Energy
|
|
|
* Laboratory.
|
|
|
*
|
|
|
* This method is a modified version of stest.c found at
|
|
|
* http://rredc.nrel.gov/solar/codes_algs/solpos/
|
|
|
*/
|
|
|
void solarzen_az(double ut, double gdlat, double glon, double * szen, double * saz)
|
|
|
{
|
|
|
struct posdata pd, *pdat; /* declare a posdata struct and a pointer for it */
|
|
|
|
|
|
int iyr = 0, imd = 0, ihm = 0, ics = 0;
|
|
|
int month = 0;
|
|
|
int day = 0;
|
|
|
int hour = 0;
|
|
|
int min = 0;
|
|
|
int sec = 0;
|
|
|
|
|
|
|
|
|
long retval; /* to capture S_solpos return codes */
|
|
|
pdat = &pd; /* point to the structure for convenience */
|
|
|
S_init (pdat);
|
|
|
|
|
|
|
|
|
/* convert UT */
|
|
|
dinvmadptr(ut, &iyr, &imd, &ihm, &ics);
|
|
|
month = imd/100;
|
|
|
day = imd - 100*month;
|
|
|
hour = ihm/100;
|
|
|
min = ihm - 100*hour;
|
|
|
sec = ics/100;
|
|
|
|
|
|
pdat->longitude = glon;
|
|
|
pdat->latitude = gdlat;
|
|
|
pdat->timezone = 0.0; /* Since we always use UT. */
|
|
|
|
|
|
|
|
|
pdat->year = iyr;
|
|
|
pdat->daynum = madGetDayno(iyr, month, day);
|
|
|
pdat->hour = hour;
|
|
|
pdat->minute = min;
|
|
|
pdat->second = sec;
|
|
|
|
|
|
retval = S_solpos (pdat); /* S_solpos function call */
|
|
|
|
|
|
/* check whether error occured */
|
|
|
if (retval != 0)
|
|
|
{
|
|
|
*szen = missing;
|
|
|
*saz = missing;
|
|
|
}
|
|
|
else
|
|
|
{
|
|
|
*szen = pdat->zenetr;
|
|
|
*saz = pdat->azim;
|
|
|
}
|
|
|
}
|
|
|
|
|
|
|
|
|
/***********************************************************************
|
|
|
*
|
|
|
* solardist calculates the distance in km from the center of the earth
|
|
|
* to the center of the sun at time ut.
|
|
|
*
|
|
|
*
|
|
|
* input parameters
|
|
|
* double ut - Universal time in seconds since 1950
|
|
|
*
|
|
|
* returns double - distance in km from the center of the earth
|
|
|
* to the center of the sun at time ut
|
|
|
*
|
|
|
* This method is taken from "Practical Astronomy with Your
|
|
|
* Calculator" 2nd edition, Peter Duffett-Smith, p. 80-87.
|
|
|
*/
|
|
|
double solardist(double ut)
|
|
|
{
|
|
|
/* constants defined in Practical Astronomy */
|
|
|
const double daysPerYear = 365.2422;
|
|
|
const double eclipLongEpic = 278.83354;
|
|
|
const double eclipLongPer = 282.596403;
|
|
|
const double eccOrbit = 0.016718;
|
|
|
const double semiMajAxis = 1.495985E8;
|
|
|
|
|
|
/* We use epic at 1950, Practical Astronomy */
|
|
|
/* uses epic at 1980, so conversion needed */
|
|
|
const double sec1950_1980 = 946684800.0;
|
|
|
|
|
|
double numDays = 0.0;
|
|
|
double numDeg = 0.0;
|
|
|
double meanAnom = 0.0;
|
|
|
double eqCenter = 0.0;
|
|
|
double trueAnom = 0.0;
|
|
|
double distance = 0.0;
|
|
|
|
|
|
/* Step one and two in book (p. 83) - computes days since 1/1/1980 */
|
|
|
numDays = (ut-sec1950_1980)/(3600.0*24.0);
|
|
|
|
|
|
/* Step three in book (p. 83) */
|
|
|
numDeg = (360.0*numDays)/daysPerYear;
|
|
|
while (numDeg < 0.0) numDeg += 360.0;
|
|
|
while (numDeg > 360.0) numDeg -= 360.0;
|
|
|
|
|
|
/* Step four in book (p. 83) */
|
|
|
meanAnom = numDeg + eclipLongEpic - eclipLongPer;
|
|
|
while (meanAnom < 0.0) meanAnom += 360.0;
|
|
|
while (meanAnom > 360.0) meanAnom -= 360.0;
|
|
|
|
|
|
/* Step five in book (p. 83) */
|
|
|
eqCenter = (360.0*eccOrbit*sin(meanAnom/57.297))/3.1415;
|
|
|
|
|
|
trueAnom = meanAnom + eqCenter;
|
|
|
while (trueAnom < 0.0) trueAnom += 360.0;
|
|
|
while (trueAnom > 360.0) trueAnom -= 360.0;
|
|
|
|
|
|
/* from section 44, p. 87 */
|
|
|
distance = semiMajAxis*(1-pow(eccOrbit, 2.0));
|
|
|
distance = distance/(1+(eccOrbit*cos(trueAnom/57.297)));
|
|
|
|
|
|
return (distance);
|
|
|
}
|
|
|
|
|
|
|
|
|
/***********************************************************************
|
|
|
*
|
|
|
* shadowheight calculates the distance directly above any gdlat and glon
|
|
|
* for a given UT in km at which the earth's shadow terminates.
|
|
|
* Will be 0.0 on dayside of earth.
|
|
|
*
|
|
|
*
|
|
|
* input parameters
|
|
|
* double ut - Universal time in seconds since 1950
|
|
|
* double gdlat - geodetic latitude in degrees
|
|
|
* double glon - geodetic longitude in degrees
|
|
|
*
|
|
|
* returns double - distance directly above any gdlat and glon
|
|
|
* for a given UT in km at which the earth's shadow terminates.
|
|
|
* Will be 0.0 on dayside of earth.
|
|
|
*
|
|
|
* This method uses the results of solarzen and solardist to create a simple
|
|
|
* cone on a sphere model of the earth's shadow. Shadow height is defined as
|
|
|
* the lowest elevation at which any part of the sun can be seen. No atmospheric
|
|
|
* bending of light is included. The radius of the earth is calculated at the
|
|
|
* tan point of the sun's rays.
|
|
|
*
|
|
|
* Algorithm:
|
|
|
*
|
|
|
* Solar Zenith Angle = Z
|
|
|
* Solar Azimuth =Az (0 = North, 90 = East)
|
|
|
*
|
|
|
* Get latitude of tangent point of sun's rays:
|
|
|
*
|
|
|
* = gdlat + cos(Az)*(Z-90.0)
|
|
|
*
|
|
|
* Get earthRadius at that point using convrt at gdlat = 0 (sea level)
|
|
|
*
|
|
|
* ConeHalfAngle = C = atan((sunRadius - earthRadius)/soldist)
|
|
|
*
|
|
|
*
|
|
|
* (cos Z tan C + 1 - sin Z)
|
|
|
* Shadowheight = earthRadius -------------------------
|
|
|
* (sin Z - cos Z tan C)
|
|
|
*
|
|
|
* Daytime (Shadowheight = 0) if numerator negitive or if
|
|
|
* Z <= 91.0.
|
|
|
*/
|
|
|
double shadowheight(double ut, double gdlat, double glon)
|
|
|
{
|
|
|
/* constants in km */
|
|
|
const double sunRadius = 695500.0;
|
|
|
|
|
|
double soldist = 0.0;
|
|
|
double szen = 0.0;
|
|
|
double saz = 0.0;
|
|
|
double coneHAng = 0.0;
|
|
|
double numer = 0.0;
|
|
|
double denom = 0.0;
|
|
|
|
|
|
/* used to get earth radius at sun tangent point */
|
|
|
double dellat = 0.0; /* change in latitude from meas point to */
|
|
|
/* sun tan point */
|
|
|
double tan_gdlat = 0.0; /* latitude of sun tangent crossing */
|
|
|
double tan_gdalt = 0.0; /* sea level sun tangent assumed */
|
|
|
double gclat = 0.0; /* don't care about geocentric lat */
|
|
|
double rkmp = 0.0; /* radius of the earth at the sun tangent */
|
|
|
|
|
|
/* get the solar zenith */
|
|
|
/* if less than 91.0, return 0.0 (dayside) */
|
|
|
solarzen_az(ut, gdlat, glon, &szen, &saz);
|
|
|
if (szen <= 91.0)
|
|
|
return (0.0);
|
|
|
|
|
|
/* find the solar tangent lat */
|
|
|
dellat =cos(saz/57.297)*(szen-90.0);
|
|
|
tan_gdlat = gdlat + dellat;
|
|
|
|
|
|
if (tan_gdlat > 90.0)
|
|
|
tan_gdlat = tan_gdlat - 2*(tan_gdlat - 90.0);
|
|
|
if (tan_gdlat < -90.0)
|
|
|
tan_gdlat = tan_gdlat - 2*(tan_gdlat + 90.0);
|
|
|
|
|
|
/* find the earth's radius at that point */
|
|
|
convrt(1, &tan_gdlat, &tan_gdalt, &gclat, &rkmp);
|
|
|
|
|
|
|
|
|
/* get distance to sun at ut */
|
|
|
soldist = solardist(ut);
|
|
|
|
|
|
/* get cone 1/2 angle */
|
|
|
coneHAng = atan((sunRadius - rkmp)/soldist);
|
|
|
|
|
|
/* get numerator */
|
|
|
numer = cos(szen/57.297)*tan(coneHAng/57.297) + 1 - sin(szen/57.297);
|
|
|
|
|
|
if (numer < 0.0) /* sun still visible on ground */
|
|
|
return (0.0);
|
|
|
|
|
|
denom = sin(szen/57.297) - cos(szen/57.297)*tan(coneHAng/57.297);
|
|
|
|
|
|
return (rkmp*numer/denom);
|
|
|
}
|
|
|
|
|
|
|
|
|
/***********************************************************************
|
|
|
*
|
|
|
* sunrise_set calculates the time UT of ionospheric sunrise
|
|
|
* and sunset.
|
|
|
*
|
|
|
*
|
|
|
* input parameters
|
|
|
* double ut - Universal time in seconds since 1950
|
|
|
* double gdlat - geodetic latitude in degrees
|
|
|
* double glon - longitude in degrees
|
|
|
* double gdalt - geodetic altitude in km
|
|
|
* double * sunrise - pointer to double allocated by user to be
|
|
|
* set to sunrise time UT
|
|
|
* double * sunset - pointer to double allocated by user to be
|
|
|
* set to sunset time UT
|
|
|
*
|
|
|
* returns void
|
|
|
*
|
|
|
* If either sunrise or sunset not found, set to missing.
|
|
|
* All times in seconds since 1/1/1950
|
|
|
*
|
|
|
* Algorithm:
|
|
|
*
|
|
|
* Depends on shadowheight calculation, so limitations discussed
|
|
|
* there apply (atmospheric bending of light ignored).
|
|
|
*
|
|
|
* If glon < 0:
|
|
|
* solar midnight = UT - glon*24/360
|
|
|
* solar noon = UT + 12 - glon*24/360
|
|
|
*
|
|
|
* If sun is not set at solar midnight (defined by shadowheight > gdalt),
|
|
|
* or sun not up at solar noon, check if any difference between 0 and 24 UT.
|
|
|
* If so, find only one of sunrise and sunset as described below, and set
|
|
|
* the other to missing. If not, return missing for both, because that point
|
|
|
* is either in the sun or in the shadow all day. The user must determine
|
|
|
* which by comparing shadowheight and gdalt. Otherwise, seek sunrise
|
|
|
* between solar midnight and solar noon, slice remaining time in half each
|
|
|
* guess. Stop when time step less than 1 minute.
|
|
|
*
|
|
|
* If sun is up at 0 UT that day: Seek sunset between 0.0 and
|
|
|
* solar midnight as above.
|
|
|
* Else: Seek sunset between solar noon and 24.0 as above.
|
|
|
*
|
|
|
* Else if glon > 0:
|
|
|
* solar midnight = UT + 24 - glon*24/360
|
|
|
* solar noon = UT + 12 - glon*24/360
|
|
|
*
|
|
|
* If sun is not set at solar midnight (defined by shadowheight > gdalt),
|
|
|
* or sun not up at solar noon, check if any difference between 0 and 24 UT.
|
|
|
* If so, find only one of sunrise and sunset as described below, and set
|
|
|
* the other to missing. If not, return missing for both, because that point
|
|
|
* is either in the sun or in the shadow all day. The user must determine
|
|
|
* which by comparing shadowheight and gdalt. Otherwise, seek sunset
|
|
|
* between solar noon and solar midnight, slice remaining time in half each
|
|
|
* guess. Stop when time step less than 1 minute.
|
|
|
*
|
|
|
* If sun is up at 0 UT that day: Seek sunrise between solar
|
|
|
* midnight and 24.0 as above.
|
|
|
*
|
|
|
* Else: Seek sunrise between 0.0 and solar noon as above.
|
|
|
*
|
|
|
* Note: day is divided 11 times to ensure greater than one minute resolution.
|
|
|
*/
|
|
|
void sunrise_set(double ut,
|
|
|
double gdlat,
|
|
|
double glon,
|
|
|
double gdalt,
|
|
|
double * sunrise,
|
|
|
double * sunset)
|
|
|
{
|
|
|
double startUT = 0.0; /* time UT at UT hour 0 */
|
|
|
double endUT = 0.0; /* time UT at UT hour 24 */
|
|
|
double solmidnight = 0.0; /* solar midnight */
|
|
|
double solnoon = 0.0; /* solar noon */
|
|
|
double startTime = 0.0; /* loop start time */
|
|
|
double endTime = 0.0; /* loop end time */
|
|
|
double tempTime = 0.0; /* loop temp time */
|
|
|
|
|
|
int iyr = 0, imd = 0, ihm = 0, ics = 0;
|
|
|
int i = 0;
|
|
|
|
|
|
int upAt0UT = 0; /* flag to indicate sun up at 0 UT */
|
|
|
int upAt24UT = 0; /* flag to indicate sun up at 24 UT */
|
|
|
|
|
|
|
|
|
/* first find UT at beginning and end of day */
|
|
|
dinvmadptr(ut, &iyr, &imd, &ihm, &ics);
|
|
|
|
|
|
startUT = dmadptr(iyr, imd, 0, 0);
|
|
|
endUT = startUT + 24.0*3600.0;
|
|
|
|
|
|
/* make sure glon between -180 and 180 */
|
|
|
while (glon > 180.0) glon -= 360.0;
|
|
|
while (glon < -180.0) glon += 360.0;
|
|
|
|
|
|
if (glon < 0.0 ) /* west of 0 degrees lat */
|
|
|
{
|
|
|
solmidnight = startUT - glon*24.0*3600.0/360.0;
|
|
|
solnoon = startUT - glon*24.0*3600.0/360.0 + 12.0*3600.0;
|
|
|
|
|
|
/* check for normal sunrise/sunset */
|
|
|
if (shadowheight(solmidnight, gdlat, glon) > gdalt &&
|
|
|
shadowheight(solnoon, gdlat, glon) == 0.0)
|
|
|
{
|
|
|
/* sunrise must occur between solmidnight and solnoon */
|
|
|
startTime = solmidnight;
|
|
|
endTime = solnoon;
|
|
|
for (i=0; i<11; i++)
|
|
|
{
|
|
|
tempTime = (startTime + endTime)/2.0;
|
|
|
if (shadowheight(tempTime, gdlat, glon) < gdalt)
|
|
|
{
|
|
|
/* sun has risen at temp time */
|
|
|
endTime = tempTime;
|
|
|
}
|
|
|
else /* sun hasn't risen yet */
|
|
|
startTime = tempTime;
|
|
|
}
|
|
|
*sunrise = (startTime + endTime)/2.0;
|
|
|
|
|
|
/* check if sun is up at 0 UT */
|
|
|
if (shadowheight(startUT, gdlat, glon) < gdalt)
|
|
|
{
|
|
|
/* seek sunset between startUT and solmidnight */
|
|
|
startTime = startUT;
|
|
|
endTime = solmidnight;
|
|
|
for (i=0; i<11; i++)
|
|
|
{
|
|
|
tempTime = (startTime + endTime)/2.0;
|
|
|
if (shadowheight(tempTime, gdlat, glon) < gdalt)
|
|
|
{
|
|
|
/* sun still up at temp time */
|
|
|
startTime = tempTime;
|
|
|
}
|
|
|
else /* sun has already set */
|
|
|
endTime = tempTime;
|
|
|
}
|
|
|
*sunset = (startTime + endTime)/2.0;
|
|
|
}
|
|
|
else
|
|
|
{
|
|
|
/* seek sunset between solnoon and endUT */
|
|
|
startTime = solnoon;
|
|
|
endTime = endUT;
|
|
|
for (i=0; i<11; i++)
|
|
|
{
|
|
|
tempTime = (startTime + endTime)/2.0;
|
|
|
if (shadowheight(tempTime, gdlat, glon) < gdalt)
|
|
|
{
|
|
|
/* sun still up at temp time */
|
|
|
startTime = tempTime;
|
|
|
}
|
|
|
else /* sun has already set */
|
|
|
endTime = tempTime;
|
|
|
}
|
|
|
*sunset = (startTime + endTime)/2.0;
|
|
|
}
|
|
|
}
|
|
|
else /* not normal - see if any transition occured */
|
|
|
{
|
|
|
upAt0UT = shadowheight(startUT, gdlat, glon) < gdalt;
|
|
|
upAt24UT = shadowheight(endUT, gdlat, glon) < gdalt;
|
|
|
|
|
|
if (upAt0UT == upAt24UT) /* no transition */
|
|
|
{
|
|
|
*sunrise = missing;
|
|
|
*sunset = missing;
|
|
|
}
|
|
|
else if (upAt0UT == 0 && upAt24UT == 1)
|
|
|
{
|
|
|
/* find sunrise only */
|
|
|
startTime = startUT;
|
|
|
endTime = endUT;
|
|
|
for (i=0; i<11; i++)
|
|
|
{
|
|
|
tempTime = (startTime + endTime)/2.0;
|
|
|
if (shadowheight(tempTime, gdlat, glon) < gdalt)
|
|
|
{
|
|
|
/* sun has risen at temp time */
|
|
|
endTime = tempTime;
|
|
|
}
|
|
|
else /* sun hasn't risen yet */
|
|
|
startTime = tempTime;
|
|
|
}
|
|
|
*sunrise = (startTime + endTime)/2.0;
|
|
|
*sunset = missing;
|
|
|
}
|
|
|
else /* sunset only */
|
|
|
{
|
|
|
startTime = startUT;
|
|
|
endTime = endUT;
|
|
|
for (i=0; i<11; i++)
|
|
|
{
|
|
|
tempTime = (startTime + endTime)/2.0;
|
|
|
if (shadowheight(tempTime, gdlat, glon) < gdalt)
|
|
|
{
|
|
|
/* sun still up at temp time */
|
|
|
startTime = tempTime;
|
|
|
}
|
|
|
else /* sun has set */
|
|
|
endTime = tempTime;
|
|
|
}
|
|
|
*sunset = (startTime + endTime)/2.0;
|
|
|
*sunrise = missing;
|
|
|
}
|
|
|
}
|
|
|
}
|
|
|
else /* west of 0 degrees lat */
|
|
|
{
|
|
|
solmidnight = startUT - glon*24.0*3600.0/360.0 + 24.0*3600.0;
|
|
|
solnoon = startUT - glon*24.0*3600.0/360.0 + 12.0*3600.0;
|
|
|
|
|
|
/* check for normal sunrise/sunset */
|
|
|
if (shadowheight(solmidnight, gdlat, glon) > gdalt &&
|
|
|
shadowheight(solnoon, gdlat, glon) == 0.0)
|
|
|
{
|
|
|
/* sunset must occur between solnoon and solmidnight */
|
|
|
startTime = solnoon;
|
|
|
endTime = solmidnight;
|
|
|
for (i=0; i<11; i++)
|
|
|
{
|
|
|
tempTime = (startTime + endTime)/2.0;
|
|
|
if (shadowheight(tempTime, gdlat, glon) < gdalt)
|
|
|
{
|
|
|
/* sun has not set yet */
|
|
|
startTime = tempTime;
|
|
|
}
|
|
|
else /* sun already */
|
|
|
endTime = tempTime;
|
|
|
}
|
|
|
*sunset = (startTime + endTime)/2.0;
|
|
|
|
|
|
/* check if sun is up at 0 UT */
|
|
|
if (shadowheight(startUT, gdlat, glon) < gdalt)
|
|
|
{
|
|
|
/* seek sunrise between solmidnight and endUT */
|
|
|
startTime = solmidnight;
|
|
|
endTime = endUT;
|
|
|
for (i=0; i<11; i++)
|
|
|
{
|
|
|
tempTime = (startTime + endTime)/2.0;
|
|
|
if (shadowheight(tempTime, gdlat, glon) < gdalt)
|
|
|
{
|
|
|
/* sun already risen */
|
|
|
endTime = tempTime;
|
|
|
}
|
|
|
else /* sun not yet risen */
|
|
|
startTime = tempTime;
|
|
|
}
|
|
|
*sunrise = (startTime + endTime)/2.0;
|
|
|
}
|
|
|
else
|
|
|
{
|
|
|
/* seek sunrise between startUT and solnoon */
|
|
|
startTime = startUT;
|
|
|
endTime = solnoon;
|
|
|
for (i=0; i<11; i++)
|
|
|
{
|
|
|
tempTime = (startTime + endTime)/2.0;
|
|
|
if (shadowheight(tempTime, gdlat, glon) < gdalt)
|
|
|
{
|
|
|
/* sun already risen */
|
|
|
endTime = tempTime;
|
|
|
}
|
|
|
else /* sun not yet risen */
|
|
|
startTime = tempTime;
|
|
|
}
|
|
|
*sunrise = (startTime + endTime)/2.0;
|
|
|
}
|
|
|
}
|
|
|
else /* not normal - see if any transition occured */
|
|
|
{
|
|
|
upAt0UT = shadowheight(startUT, gdlat, glon) < gdalt;
|
|
|
upAt24UT = shadowheight(endUT, gdlat, glon) < gdalt;
|
|
|
|
|
|
if (upAt0UT == upAt24UT) /* no transition */
|
|
|
{
|
|
|
*sunrise = missing;
|
|
|
*sunset = missing;
|
|
|
}
|
|
|
else if (upAt0UT == 0 && upAt24UT == 1)
|
|
|
{
|
|
|
/* find sunrise only */
|
|
|
startTime = startUT;
|
|
|
endTime = endUT;
|
|
|
for (i=0; i<11; i++)
|
|
|
{
|
|
|
tempTime = (startTime + endTime)/2.0;
|
|
|
if (shadowheight(tempTime, gdlat, glon) < gdalt)
|
|
|
{
|
|
|
/* sun has risen at temp time */
|
|
|
endTime = tempTime;
|
|
|
}
|
|
|
else /* sun hasn't risen yet */
|
|
|
startTime = tempTime;
|
|
|
}
|
|
|
*sunrise = (startTime + endTime)/2.0;
|
|
|
*sunset = missing;
|
|
|
}
|
|
|
else /* sunset only */
|
|
|
{
|
|
|
startTime = startUT;
|
|
|
endTime = endUT;
|
|
|
for (i=0; i<11; i++)
|
|
|
{
|
|
|
tempTime = (startTime + endTime)/2.0;
|
|
|
if (shadowheight(tempTime, gdlat, glon) < gdalt)
|
|
|
{
|
|
|
/* sun still up at temp time */
|
|
|
startTime = tempTime;
|
|
|
}
|
|
|
else /* sun has set */
|
|
|
endTime = tempTime;
|
|
|
}
|
|
|
*sunset = (startTime + endTime)/2.0;
|
|
|
*sunrise = missing;
|
|
|
}
|
|
|
|
|
|
}/* end check for normal sunrise/sunset */
|
|
|
|
|
|
} /* end if glon < 0 */
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
