Index: /trunk/Ohana/src/skycalc/Makefile =================================================================== --- /trunk/Ohana/src/skycalc/Makefile (revision 2496) +++ /trunk/Ohana/src/skycalc/Makefile (revision 2496) @@ -0,0 +1,71 @@ +include ../../Configure +HOME = . + +SRC = $(HOME)/src +LIB = $(HOME)/lib +INC = $(HOME)/include +MAN = $(HOME)/doc +BIN = $(HOME)/bin +DESTBIN = $(LBIN) +DESTLIB = $(LLIB) +DESTINC = $(LINC) +DESTMAN = $(LMAN) + +INCS = -I$(INC) -I$(LINC) -I$(LINC)/$(ARCH) +LIBS = -L$(LLIB) +CFLAGS = -g -o $*.$(ARCH).o $(INCS) +CCFLAGS = $(INCS) $(LIBS) + +default: skycalc + +skycalc: $(DESTLIB)/libskycalc.a + +install: $(DESTINC)/skycalc.h $(DESTLIB)/libskycalc.a + +OBJ = $(SRC)/time.$(ARCH).o $(SRC)/geometry.$(ARCH).o \ + $(SRC)/astro.$(ARCH).o $(SRC)/sun.$(ARCH).o \ + $(SRC)/moon.$(ARCH).o $(SRC)/misc.$(ARCH).o + +$(OBJ): $(INC)/skycalc.h $(INC)/skycalc_internal.h + +$(LIB)/libskycalc.$(ARCH).a: $(OBJ) + rm -f $(LIB)/libskycalc.$(ARCH).a + ar rcv $(LIB)/libskycalc.$(ARCH).a $(OBJ) + ranlib $(LIB)/libskycalc.$(ARCH).a + +$(DESTLIB)/libskycalc.a: $(LIB)/libskycalc.$(ARCH).a + @if [ ! -d $(DESTLIB) ]; then mkdir -p $(DESTLIB); fi + rm -f $(DESTLIB)/libskycalc.a + cp $(LIB)/libskycalc.$(ARCH).a $(DESTLIB)/libskycalc.a + +$(DESTINC)/skycalc.h: $(INC)/skycalc.h + @if [ ! -d $(DESTINC) ]; then mkdir -p $(DESTINC); fi + rm -f $(DESTINC)/skycalc.h + cp $(INC)/skycalc.h $(DESTINC)/ + +test: + $(CC) $(CCFLAGS) -g -o skylib $(SRC)/skylib.c -lskycalc -lm + +dusktime: skycalc + $(CC) $(CCFLAGS) -g -o $(BIN)/dusktime.$(ARCH) $(SRC)/dusktime.c -lskycalc -lm -lohana + cp $(BIN)/dusktime.$(ARCH) $(DESTBIN)/dusktime + +moondata: skycalc + $(CC) $(CCFLAGS) -g -o $(BIN)/moondata.$(ARCH) $(SRC)/moondata.c -lskycalc -lm -lohana + cp $(BIN)/moondata.$(ARCH) $(DESTBIN)/moondata + +sundata: skycalc + $(CC) $(CCFLAGS) -g -o $(BIN)/sundata.$(ARCH) $(SRC)/sundata.c -lskycalc -lm -lohana + cp $(BIN)/sundata.$(ARCH) $(DESTBIN)/sundata + +clean: + rm -f `find . -name "*.o"` + rm -f `find . -name "*.a"` + rm -f `find . -name "*~"` + rm -f `find . -name "#*"` + +.SUFFIXES: .$(ARCH).o + +.c.$(ARCH).o: + $(CC) $(CFLAGS) -c $< + Index: /trunk/Ohana/src/skycalc/doc/skycalc.txt =================================================================== --- /trunk/Ohana/src/skycalc/doc/skycalc.txt (revision 2496) +++ /trunk/Ohana/src/skycalc/doc/skycalc.txt (revision 2496) @@ -0,0 +1,92 @@ + +functions in time.c: + +int get_sys_date (struct SC_date_time *date) + set SC date & time from system clock + +double date_to_jd (struct SC_date_time date) + convert JD to SC date & time + +void jd_to_date (double jdin, struct SC_date_time *date) + convert SC date & time to JD + +double lst (double jd, double longit) + convert JD & longitude to LST + +double adj_time (double x) + force time domain to be -12h and 12h + + +functions in geometry.c: + +void xyz_cel (double x, double y, double z, double *r, double *d) + converts a coordinate triplet back to a standard ra and dec + +double atan_circ (double x, double y) + returns radian angle 0 to 2pi for coords x, y -- get that quadrant right !! */ + +double altit (double dec, double ha, double lat, double *az) + returns altitude (degrees) for dec, ha, lat (degree) + also computes and returns azimuth through pointer argument. + +double ha_alt (double dec, double lat, double alt) + returns hour angle at which object at dec is at altitude alt. + If object is never at this altitude, signals with special + return values 1000 (always higher) and -1000 (always lower). + +void min_max_alt (double lat, double dec, double *min, double *max) + computes minimum and maximum altitude for a given dec and latitude. + +double circulo (double x) + force domain to be 0 - 360 degrees + + +functions in astro.c: + +void precrot (double rorig, double dorig, double orig_epoch, double final_epoch, double *rf, double *df) { + accurate precession method from rorig, dorig, orig_epoch + +void geocent (double geolong, double geolat, double height, double x_geo, double y_geo, double z_geo) + return geocentric coordinate from geodetic data + +void eclrot(double jd, double x, double y, double z) + rotates ecliptic rectangular coords x, y, z to equatorial + +double etcorr (double jd) + return correction to jd for shifts in ET or TDT + +void set_zenith (struct SC_date_time date, double lat, double longit, double epoch, double *ra, double *dec) + return zenith coords for given date, time, epoch, position + + +functions is sun.c: + +void lpsun (double jd, double ra, double dec) + return ra & dec for sun on given JD (low precision) + +double jd_sun_alt (double alt, double jdguess, double lat, double longit) + return JD for sun at given altitude, start with jdguess + +double sunset_tonight (struct SC_date_time date, double lat, double longit, double elev) + return JD for sunset before closest midnight to date & time + +double sunrise_tonight (struct SC_date_time date, double lat, double longit, double elev) + return JD for sunrise after closest midnight to date & time + + +functions in moon.c: + +void accumoon (double jd, double geolat, double lst, double elevsea, + double *geora, double *geodec, double *geodist, + double *topora, double *topodec, double *topodist) + calculate position of moon given jd, lat, etc. + +double jd_moon_alt (double alt, double jdguess, double lat, double longit, double elevsea) + return JD for moon at given altitude, start with jdguess + +double moonset_tonight (struct SC_date_time date, double lat, double longit, double elevsea, double elev) { + Given site position, return Moonset for closest midnight + +double moonrise_tonight (struct SC_date_time date, double lat, double longit, double elevsea, double elev) { + Given site position, return Moonrise for closest midnight + Index: /trunk/Ohana/src/skycalc/include/skycalc.h =================================================================== --- /trunk/Ohana/src/skycalc/include/skycalc.h (revision 2496) +++ /trunk/Ohana/src/skycalc/include/skycalc.h (revision 2496) @@ -0,0 +1,23 @@ + +/* header for skycalc library function calls */ + +struct SC_coord { + short sign; /* carry sign explicitly since -0 not neg. */ + double hh; + double mm; + double ss; +}; + +struct SC_date_time { + short y; + short mo; + short d; + short h; + short mn; + float s; +}; + +double sunset_tonight (struct SC_date_time date, double lat, double longit, double elev); +double sunrise_tonight (struct SC_date_time date, double lat, double longit, double elev); +double moonset_tonight (struct SC_date_time date, double lat, double longit, double elevsea, double elev); +double moonrise_tonight (struct SC_date_time date, double lat, double longit, double elevsea, double elev); Index: /trunk/Ohana/src/skycalc/include/skycalc_internal.h =================================================================== --- /trunk/Ohana/src/skycalc/include/skycalc_internal.h (revision 2496) +++ /trunk/Ohana/src/skycalc/include/skycalc_internal.h (revision 2496) @@ -0,0 +1,91 @@ + +/* header for use by the skycalc library files, not needed for external calls */ + +# include +# include +# include +# include +# include +# include +# include + +/* a couple of the system-dependent magic numbers are defined here */ + +#define SYS_CLOCK_OK 1 /* 1 means ANSI-standard time libraries do work, + 2 means they don't. This is used by compiler switches in file 5 and + the main program. */ + +#define LOG_FILES_OK 1 /* 1 means that log files are enabled. + Any other value means they're not. */ + +#define MAX_OBJECTS 500 +#define MINSHORT -32767 /* min, max short integers and double precision */ +#define MAXSHORT 32767 +#define MAXDOUBLE 1.0e38 +#define MINDOUBLE -1.0e38 +#define BUFSIZE 150 + +/* some (not all) physical, mathematical, and astronomical constants + used are defined here. */ + +#define PI 3.14159265358979 +#define ARCSEC_IN_RADIAN 206264.8062471 +#define DEG_IN_RADIAN 57.2957795130823 +#define HRS_IN_RADIAN 3.819718634205 +#define KMS_AUDAY 1731.45683633 /* km per sec in 1 AU/day */ +#define SS_MASS 1.00134198 /* solar system mass in solar units */ +#define J2000 2451545. /* Julian date at standard epoch */ +#define SEC_IN_DAY 86400. +#define FLATTEN 0.003352813 /* flattening of earth, 1/298.257 */ +#define EQUAT_RAD 6378137. /* equatorial radius of earth, meters */ +#define ASTRO_UNIT 1.4959787066e11 /* 1 AU in meters */ +#define RSUN 6.96000e8 /* IAU 1976 recom. solar radius, meters */ +#define RMOON 1.738e6 /* IAU 1976 recom. lunar radius, meters */ +#define PLANET_TOL 3. /* flag if nearer than 3 degrees + to a major planet ... */ +#define KZEN 0.172 /* zenith extinction, mag, for use + in lunar sky brightness calculations. */ +#define FIRSTJD 2415387. /* 1901 Jan 1 -- calendrical limit */ +#define LASTJD 2488070. /* 2099 Dec 31 */ + +/* MAGIC NUMBERS which might depend on how accurately double- + precision floating point is handled on your machine ... */ + +#define EARTH_DIFF 0.05 /* used in numerical + differentiation to find earth velocity -- this value gives + about 8 digits of numerical accuracy on the VAX, but is + about 3 orders of magnitude larger than the value where roundoff + errors become apparent. */ + +#define MIDN_TOL 0.00001 /* this is no longer + used -- it was formerly + how close (in days) a julian date has to be to midnight + before a warning flag is printed for the reader. VAX + double precision renders a Julian date considerably + more accurately than this. The day and date are now based + on the same rounding of the julian date, so they should + always agree. */ + +/** prototypes **/ +int get_sys_date (struct SC_date_time *date); +double date_to_jd (struct SC_date_time date); +void jd_to_date (double jdin, struct SC_date_time *date); +double lst (double jd, double longit); +double adj_time (double x); +void xyz_cel (double x, double y, double z, double *r, double *d); +double atan_circ (double x, double y); +double altit (double dec, double ha, double lat, double *az); +double ha_alt (double dec, double lat, double alt); +void min_max_alt (double lat, double dec, double *min, double *max); +double circulo (double x); +void precrot (double rorig, double dorig, double orig_epoch, double final_epoch, double *rf, double *df); +void geocent (double geolong, double geolat, double height, double *x_geo, double *y_geo, double *z_geo); +void eclrot(double jd, double *x, double *y, double *z); +double etcorr (double jd); +void set_zenith (struct SC_date_time date, double lat, double longit, double epoch, double *ra, double *dec); +void lpsun (double jd, double *ra, double *dec); +double jd_sun_alt (double alt, double jdguess, double lat, double longit); +void accumoon (double jd, double geolat, double lst, double elevsea, + double *geora, double *geodec, double *geodist, + double *topora, double *topodec, double *topodist); +double jd_moon_alt (double alt, double jdguess, double lat, double longit, double elevsea); Index: /trunk/Ohana/src/skycalc/src/astro.c =================================================================== --- /trunk/Ohana/src/skycalc/src/astro.c (revision 2496) +++ /trunk/Ohana/src/skycalc/src/astro.c (revision 2496) @@ -0,0 +1,192 @@ +# include + +/* Takes a coordinate pair and precesses it using matrix procedures + as outlined in Taff's Computational Spherical Astronomy book. + This is the so-called 'rigorous' method which should give very + accurate answers all over the sky over an interval of several + centuries. Naked eye accuracy holds to ancient times, too. + Precession constants used are the new IAU1976 -- the 'J2000' + system. Angles in degrees, epochs in years */ + +void precrot (double rorig, double dorig, double orig_epoch, double final_epoch, double *rf, double *df) { + + double ti, tf, zeta, z, theta; /* all as per Taff */ + double cosz, coszeta, costheta, sinz, sinzeta, sintheta; /* ftns */ + double p11, p12, p13, p21, p22, p23, p31, p32, p33; + /* elements of the rotation matrix */ + double radian_ra, radian_dec; + double orig_x, orig_y, orig_z; + double fin_x, fin_y, fin_z; /* original and final unit ectors */ + + ti = (orig_epoch - 2000.) / 100.; + tf = (final_epoch - 2000. - 100. * ti) / 100.; + + zeta = (2306.2181 + 1.39656 * ti + 0.000139 * ti * ti) * tf + + (0.30188 - 0.000344 * ti) * tf * tf + 0.017998 * tf * tf * tf; + z = zeta + (0.79280 + 0.000410 * ti) * tf * tf + 0.000205 * tf * tf * tf; + theta = (2004.3109 - 0.8533 * ti - 0.000217 * ti * ti) * tf + - (0.42665 + 0.000217 * ti) * tf * tf - 0.041833 * tf * tf * tf; + + /* convert to radians */ + + zeta = zeta / ARCSEC_IN_RADIAN; + z = z / ARCSEC_IN_RADIAN; + theta = theta / ARCSEC_IN_RADIAN; + + /* compute the necessary trig functions for speed and simplicity */ + + cosz = cos(z); + coszeta = cos(zeta); + costheta = cos(theta); + sinz = sin(z); + sinzeta = sin(zeta); + sintheta = sin(theta); + + /* compute the elements of the precession matrix */ + + p11 = coszeta * cosz * costheta - sinzeta * sinz; + p12 = -1. * sinzeta * cosz * costheta - coszeta * sinz; + p13 = -1. * cosz * sintheta; + + p21 = coszeta * sinz * costheta + sinzeta * cosz; + p22 = -1. * sinzeta * sinz * costheta + coszeta * cosz; + p23 = -1. * sinz * sintheta; + + p31 = coszeta * sintheta; + p32 = -1. * sinzeta * sintheta; + p33 = costheta; + + /* transform original coordinates */ + + radian_ra = rorig / HRS_IN_RADIAN; + radian_dec = dorig / DEG_IN_RADIAN; + + orig_x = cos(radian_dec) * cos(radian_ra); + orig_y = cos(radian_dec) *sin(radian_ra); + orig_z = sin(radian_dec); + /* (hard coded matrix multiplication ...) */ + fin_x = p11 * orig_x + p12 * orig_y + p13 * orig_z; + fin_y = p21 * orig_x + p22 * orig_y + p23 * orig_z; + fin_z = p31 * orig_x + p32 * orig_y + p33 * orig_z; + + /* convert back to spherical polar coords */ + + xyz_cel(fin_x, fin_y, fin_z, rf, df); + +} + + +/* computes the geocentric coordinates from the geodetic + (standard map-type) longitude, latitude, and height. + These are assumed to be in decimal hours, decimal degrees, and + meters respectively. Notation generally follows 1992 Astr Almanac, + p. K11 */ + +void geocent (double geolong, double geolat, double height, double *x_geo, double *y_geo, double *z_geo) { + + double denom, C_geo, S_geo; + + geolat = geolat / DEG_IN_RADIAN; + geolong = geolong / HRS_IN_RADIAN; + denom = (1. - FLATTEN) * sin(geolat); + denom = cos(geolat) * cos(geolat) + denom*denom; + C_geo = 1. / sqrt(denom); + S_geo = (1. - FLATTEN) * (1. - FLATTEN) * C_geo; + C_geo = C_geo + height / EQUAT_RAD; /* deviation from almanac + notation -- include height here. */ + S_geo = S_geo + height / EQUAT_RAD; + *x_geo = C_geo * cos(geolat) * cos(geolong); + *y_geo = C_geo * cos(geolat) * sin(geolong); + *z_geo = S_geo * sin(geolat); +} + +/* rotates ecliptic rectangular coords x, y, z to + equatorial (all assumed of date.) */ +void eclrot(double jd, double *x, double *y, double *z) { + + double incl; + double xpr,ypr,zpr; + double T; + + T = (jd - J2000) / 36525; /* centuries since J2000 */ + + incl = (23.439291 + T * (-0.0130042 - 0.00000016 * T))/DEG_IN_RADIAN; + /* 1992 Astron Almanac, p. B18, dropping the + cubic term, which is 2 milli-arcsec! */ + ypr = cos(incl) * *y - sin(incl) * *z; + zpr = sin(incl) * *y + cos(incl) * *z; + *y = ypr; + *z = zpr; + /* x remains the same. */ +} + +/* Given a julian date in 1900-2100, returns the correction + delta t which is: + TDT - UT (after 1983 and before 1993) + ET - UT (before 1983) + an extrapolated guess (after 1993). + + For dates in the past (<= 1993) the value is linearly + interpolated on 5-year intervals; for dates after the present, + an extrapolation is used, because the true value of delta t + cannot be predicted precisely. Note that TDT is essentially the + modern version of ephemeris time with a slightly cleaner + definition. + + Where the algorithm shifts there is an approximately 0.1 second + discontinuity. Also, the 5-year linear interpolation scheme can + lead to errors as large as 0.5 seconds in some cases, though + usually rather smaller. */ + +double etcorr (double jd) { + + double jd1900 = 2415019.5; + double dates[20] = {1900,1905,1910,1915,1920,1925,1930,1935,1940,1945, + 1950,1955,1960,1965,1970,1975,1980,1985,1990,1993}; + double delts[20]={-2.72,3.86,10.46,17.20,21.16,23.62,24.02,23.93,24.33,26.77, + 29.15,31.07,33.15,35.73,40.18,45.48,50.54,54.34,56.86,59.12}; + double year, delt; + short i; + + year = 1900. + (jd - 2415019.5) / 365.25; + + if(year < 1993.0 && year >= 1900.) { + i = (year - 1900) / 5; + delt = delts[i] + + ((delts[i+1] - delts[i])/(dates[i+1] - dates[i])) * (year - dates[i]); + } + + else if (year > 1993. && year < 2100.) + delt = 33.15 + (2.164e-3) * (jd - 2436935.4); /* rough extrapolation */ + + else if (year < 1900) { + delt = 0.; + } + + else if (year >= 2100.) { + delt = 180.; /* who knows? */ + } + + return (delt); + +} + +/* sets RA and DEC at zenith as defined by given time and date */ +void set_zenith (struct SC_date_time date, double lat, double longit, double epoch, double *ra, double *dec) +{ + double jd, current_epoch; + + jd = date_to_jd (date); + + if (jd < 0.) return; /* nonexistent time. */ + + *ra = lst (jd, longit); + + *dec = lat; + + current_epoch = 2000. + (jd - J2000) / 365.25; + + precrot (*ra,*dec,current_epoch,epoch,ra,dec); + +} + Index: /trunk/Ohana/src/skycalc/src/dusktime.c =================================================================== --- /trunk/Ohana/src/skycalc/src/dusktime.c (revision 2496) +++ /trunk/Ohana/src/skycalc/src/dusktime.c (revision 2496) @@ -0,0 +1,69 @@ +# include +# include +# include + +void set_site (double *longit, double *lat, double *elevsea, double *elev) { + + *longit = 10.36478; /* W longitude in decimal hours */ + *lat = 19.8267; /* N latitude in decimal degrees */ + *elevsea = 4215.; /* elevation above sea level (for absolute location) */ + *elev = 4215.; /* observatory elevation above horizon, meters */ + +} + +main (int argc, char **argv) { + + struct SC_date_time date, tmpdate; + double longit, lat, elevsea, elev; + double jdnow, jdset, jdrise, d1, d2; + unsigned long int tzero; + struct tm *stm; + + if (argc != 2) { + fprintf (stderr, "USAGE: dusktime (date)\n"); + exit (1); + } + + if (!str_to_time (argv[1], &tzero)) { + fprintf (stderr, "syntax error\n"); + exit (1); + } + stm = gmtime (&tzero); + date.y = (short) (stm->tm_year + 1900); + date.mo = (short) (stm->tm_mon + 1); + date.d = (short) (stm->tm_mday); + date.h = (short) (stm->tm_hour); + date.mn = (short) (stm->tm_min); + date.s = (float) (stm->tm_sec); + + /* + fprintf (stderr, "%4d/%02d/%02d %02d:%02d:%02f\n", date.y, date.mo, date.d, date.h, date.mn, date.s); + */ + + set_site (&longit, &lat, &elevsea, &elev); + + jdnow = date_to_jd (date); + jdset = sunset_tonight (date, lat, longit, elev); + jdrise = sunrise_tonight (date, lat, longit, elev); + + + d1 = fabs (jdnow - jdset); + d2 = fabs (jdnow - jdrise); + + if (d1 < d2) { + fprintf (stdout, "set %f %f\n", jdset, 24*60*(jdnow - jdset)); + jdnow = jdset; + } else { + fprintf (stdout, "rise %f %f\n", jdrise, 24*60*(jdrise - jdnow)); + jdnow = jdrise; + } + + /* + jd_to_date (jdnow, &tmpdate); + fprintf (stderr, "%4d/%02d/%02d %02d:%02d:%02f\n", tmpdate.y, tmpdate.mo, tmpdate.d, tmpdate.h, tmpdate.mn, tmpdate.s); + */ + exit (0); + +} + + /* set_zenith (date, lat, longit, objepoch, &objra, &objdec); */ Index: /trunk/Ohana/src/skycalc/src/geometry.c =================================================================== --- /trunk/Ohana/src/skycalc/src/geometry.c (revision 2496) +++ /trunk/Ohana/src/skycalc/src/geometry.c (revision 2496) @@ -0,0 +1,135 @@ +# include + +/* converts a coordinate triplet back to a standard ra and dec */ +void xyz_cel (double x, double y, double z, double *r, double *d) { + + /* converts a coordinate triplet back to a standard ra and dec */ + + double mod; /* modulus */ + double xy; /* component in xy plane */ + short sign; /* for determining quadrant */ + double radian_ra, radian_dec; + + /* this taken directly from pl1 routine - no acos or asin available there, + as it is in c. Easier just to copy, though */ + + mod = sqrt(x*x + y*y + z*z); + x = x / mod; + y = y / mod; + z = z / mod; /* normalize 'em explicitly first. */ + + xy = sqrt(x*x + y*y); + + if(xy < 1.0e-10) { + radian_ra = 0.; /* too close to pole */ + radian_dec = PI / 2.; + if(z < 0.) radian_dec = radian_dec * -1.; + } + else { + if(fabs(z/xy) < 3.) radian_dec = atan(z / xy); + else if (z >= 0.) radian_dec = PI / 2. - atan(xy / z); + else radian_dec = -1. * PI / 2. - atan(xy / z); + if(fabs(x) > 1.0e-10) { + if(fabs(y / x) < 3.) radian_ra = atan(y / x); + else if ((x * y ) >= 0.) radian_ra = PI / 2. - atan(x/y); + else radian_ra = -1. * PI / 2. - atan(x / y); + } + else { + radian_ra = PI / 2.; + if((x * y)<= 0.) radian_ra = radian_ra * -1.; + } + if(x <0.) radian_ra = radian_ra + PI ; + if(radian_ra < 0.) radian_ra = radian_ra + 2. * PI ; + } + + *r = radian_ra * HRS_IN_RADIAN; + *d = radian_dec * DEG_IN_RADIAN; + +} + +/* returns radian angle 0 to 2pi for coords x, y -- get that quadrant right !! */ +double atan_circ (double x, double y) { + + double theta; + + if(x == 0.) { + if(y > 0.) theta = PI / 2.; + else if(y < 0.) theta = 3.* PI / 2.; + else theta = 0.; /* x and y zero */ + } + else theta = atan(y/x); + if(x < 0.) theta = theta + PI; + if(theta < 0.) theta = theta + 2.* PI; + return(theta); +} + +/* returns altitude(degr) for dec, ha, lat (decimal degr, hr, degr); + also computes and returns azimuth through pointer argument. */ +double altit (double dec, double ha, double lat, double *az) { + + double x,y,z; + + dec = dec / DEG_IN_RADIAN; + ha = ha / HRS_IN_RADIAN; + lat = lat / DEG_IN_RADIAN; + x = DEG_IN_RADIAN * asin(cos(dec)*cos(ha)*cos(lat) + sin(dec)*sin(lat)); + y = sin(dec)*cos(lat) - cos(dec)*cos(ha)*sin(lat); /* due N comp. */ + z = -1. * cos(dec)*sin(ha); /* due east comp. */ + *az = atan_circ(y,z) * DEG_IN_RADIAN; + return(x); +} + +/* returns hour angle at which object at dec is at altitude alt. + If object is never at this altitude, signals with special + return values 1000 (always higher) and -1000 (always lower). */ +double ha_alt (double dec, double lat, double alt) { + + double x,coalt,min,max; + + min_max_alt(lat,dec,&min,&max); + if(alt < min) + return(1000.); /* flag value - always higher than asked */ + if(alt > max) + return(-1000.); /* flag for object always lower than asked */ + dec = (0.5*PI) - dec / DEG_IN_RADIAN; + lat = (0.5*PI) - lat / DEG_IN_RADIAN; + coalt = (0.5*PI) - alt / DEG_IN_RADIAN; + x = (cos(coalt) - cos(dec)*cos(lat)) / (sin(dec)*sin(lat)); + if (fabs(x) <= 1.) return(acos(x) * HRS_IN_RADIAN); + else { + printf ("Error in ha_alt ... acos(>1).\n"); + return (1000.); + } +} + +/* computes minimum and maximum altitude for a given dec and latitude. */ +void min_max_alt (double lat, double dec, double *min, double *max) { + + double x; + + lat = lat / DEG_IN_RADIAN; /* pass by value! */ + dec = dec / DEG_IN_RADIAN; + x = cos(dec)*cos(lat) + sin(dec)*sin(lat); + if (fabs(x) <= 1.) { + *max = asin(x) * DEG_IN_RADIAN; + } + else printf ("Error in min_max_alt -- arcsin(>1)\n"); + + x = sin(dec)*sin(lat) - cos(dec)*cos(lat); + if(fabs(x) <= 1.) { + *min = asin(x) * DEG_IN_RADIAN; + } + else printf ("Error in min_max_alt -- arcsin(>1)\n"); +} + +/* force domain to be 0 - 360 degrees */ +double circulo (double x) { + + /* fails for negative angles! */ + + int n; + + n = (int)(x / 360.); + return(x - 360. * n); +} + Index: /trunk/Ohana/src/skycalc/src/misc.c =================================================================== --- /trunk/Ohana/src/skycalc/src/misc.c (revision 2496) +++ /trunk/Ohana/src/skycalc/src/misc.c (revision 2496) @@ -0,0 +1,321 @@ +# include + +int dms_to_ddd (double *Value, char *string); +int str_to_radec (double *ra, double *dec, char *str1, char *str2); +int chk_time (char *line); +double sec_to_jd (unsigned long second); +unsigned long int jd_to_sec (double jd); +char *sec_to_date (unsigned long second); +unsigned long date_to_sec (char *date); + +/***** convert [-]00:00:00 to 0.0000 ****/ +int dms_to_ddd (double *Value, char *string) { + + int valid, neg, status; + double tmp, value; + char *p1, *p2, *px; + + valid = FALSE; + neg = FALSE; + stripwhite (string); + p1 = string; + px = string + strlen(string); + + if (string[0] == '-') { + valid = TRUE; + neg = TRUE; + p1 = &string[1]; + } + if (string[0] == '+') { + valid = TRUE; + neg = FALSE; + p1 = &string[1]; + } + if (isdigit(string[0])) { + valid = TRUE; + p1 = &string[0]; + } + if (!valid) { return (FALSE); } + + status = 1; + tmp = strtod (p1, &p2); + if (p2 == p1) return (FALSE); /* entry not a number: +fred */ + value = tmp; + if (p2 == px) goto escape; /* entry only number: +1.0 */ + p1 = p2 + 1; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; /* entry not a number: +1:fred */ + status = 2; + value += tmp / 60.0; + if (p2 == px) goto escape; /* entry only number: +1:1 */ + p1 = p2 + 1; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; /* entry not a number: +1:1:fred */ + value += tmp / 3600.0; + + escape: + if (neg) { + value *= -1; + } + *Value = value; + + return (status); +} + +/**********/ +int str_to_radec (double *ra, double *dec, char *str1, char *str2) { + + double Ra, Dec; + + *ra = *dec = 0; + switch (dms_to_ddd (&Ra, str1)) { + case 0: + fprintf (stderr, "syntax error in RA\n"); + return (FALSE); + case 1: + break; + case 2: + Ra = Ra * 15; + break; + } + switch (dms_to_ddd (&Dec, str2)) { + case 0: + fprintf (stderr, "syntax error in DEC\n"); + return (FALSE); + case 1: + case 2: + break; + } + *ra = Ra; + *dec = Dec; + return (TRUE); +} + +# define FORMAT_DAYS 1 +# define FORMAT_HOURS 2 +# define FORMAT_MINUTES 3 +# define FORMAT_SECONDS 4 +# define FORMAT_JD 5 +# define FORMAT_DATE 6 + +/**********/ +int chk_time (char *line) { + + char *p1, *p2; + double tmp; + int mode; + + p1 = line; + tmp = strtod (p1, &p2); + if (p2 == p1 + strlen (p1) - 1) { + if (*p2 == 'd') { + mode = FORMAT_DAYS; + } + if (*p2 == 'h') { + mode = FORMAT_HOURS; + } + if (*p2 == 'm') { + mode = FORMAT_MINUTES; + } + if (*p2 == 's') { + mode = FORMAT_SECONDS; + } + if (*p2 == 'j') { + mode = FORMAT_JD; + } + } else { + mode = FORMAT_DATE; + } + return (mode); +} + +/**********/ +int str_to_time (char *line, unsigned int *second) { + + struct timeval now; + double jd; + + if (!strcasecmp (line, "NOW")) { + gettimeofday (&now, (struct timezone *) NULL); + *second = now.tv_sec; + return (TRUE); + } + + if (!strcasecmp (line, "TODAY")) { + gettimeofday (&now, (struct timezone *) NULL); + *second = 86400 * ((int)(now.tv_sec / 86400)); + return (TRUE); + } + + switch (chk_time (line)) { + case 0: + return (FALSE); + case FORMAT_DAYS: + *second = strtod (line, 0) * 86400.0; + return (TRUE); + case FORMAT_HOURS: + *second = strtod (line, 0) * 3600.0; + return (TRUE); + case FORMAT_MINUTES: + *second = strtod (line, 0) * 60.0; + return (TRUE); + case FORMAT_SECONDS: + *second = strtod (line, 0); + return (TRUE); + case FORMAT_JD: + jd = strtod (line, 0); + *second = jd_to_sec (jd); + return (TRUE); + case FORMAT_DATE: + *second = date_to_sec (line); + return (TRUE); + } + return (FALSE); +} + +/**********/ +int str_to_dtime (char *line, double *second) { + + switch (chk_time (line)) { + case 0: + case FORMAT_JD: + case FORMAT_DATE: + return (FALSE); + case FORMAT_DAYS: + *second = strtod (line, 0) * 86400.0; + return (TRUE); + case FORMAT_HOURS: + *second = strtod (line, 0) * 3600.0; + return (TRUE); + case FORMAT_MINUTES: + *second = strtod (line, 0) * 60.0; + return (TRUE); + case FORMAT_SECONDS: + *second = strtod (line, 0); + return (TRUE); + } + return (FALSE); +} + +/**********/ +double sec_to_jd (unsigned long second) { + + double jd; + + jd = second/86400.0 + 2440587.5; + return (jd); +} + +/**********/ +unsigned long int jd_to_sec (double jd) { + + unsigned long int second; + + second = (jd - 2440587.5)*86400; + return (second); +} + +/**********/ +char *sec_to_date (unsigned long second) { + + struct tm *gmt; + char *line; + + ALLOCATE (line, char, 64); + gmt = gmtime (&second); + sprintf (line, "%4d/%02d/%02d,%02d:%02d:%02d", 1900+gmt[0].tm_year, gmt[0].tm_mon+1, gmt[0].tm_mday, gmt[0].tm_hour, gmt[0].tm_min, gmt[0].tm_sec); + return (line); + +} + +/***** date in format yyyy/mm/dd,hh:mm:ss *****/ +unsigned long date_to_sec (char *date) { + + unsigned long second; + double tmp, jd; + struct tm now; + char *p1, *p2, *px; + + p1 = date; + px = date + strlen(date); + bzero (&now, sizeof(now)); + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; + now.tm_year = tmp; + if (now.tm_year > 1000) now.tm_year -= 1900; + if (now.tm_year < 50) now.tm_year += 100; + if (p2 == px) goto escape; + p1 = p2 + 1; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; + now.tm_mon = tmp - 1; /* mon runs from 0 - 11 */ + if (p2 == px) goto escape; + p1 = p2 + 1; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; + now.tm_mday = tmp; + if (p2 == px) goto escape; + p1 = p2 + 1; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; + p1 = p2 + 1; + now.tm_hour = tmp; + if (p2 == px) goto escape; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; + now.tm_min = tmp; + if (p2 == px) goto escape; + p1 = p2 + 1; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; + now.tm_sec = tmp; + if (p2 == px) goto escape; + p1 = p2 + 1; + + escape: + jd = now.tm_mday - 32075 + (int)(1461*(1900 + now.tm_year + 4800 + (int)(((now.tm_mon+1)-14)/12))/4) + + (int)(367*((now.tm_mon+1) - 2 - (int)(((now.tm_mon+1) - 14)/12)*12)/12) + - (int)(3*(int)((1900 + now.tm_year + 4900 + (int)(((now.tm_mon+1) - 14)/12))/100)/4) - 0.5; + + second = (jd - 2440587.5)*86400 + 3600.0*now.tm_hour + now.tm_min*60.0 + now.tm_sec; + + return (second); +} + + +/* times may be in forms as: + 20040200450s (N seconds since 1970.0) + 2440900.232j (julian date) + 99/02/23,03:22:18 (date string) + (separators may be anything except space, +, -) + 99:02:15:12:23:30 + 99:02:15:12h23m30s + */ + + +/* fseek with timeout - 0.5 sec */ +int Fseek (FILE *f, long offset, int whence) { + + int status, k; + + status = fseek (f, offset, whence); + if (status == -1) { + int k; + /* fprintf (stderr, "problem seeking position: %d\n", errno); */ + for (k = 0; (k < 10) && ((status = fseek (f, 0, SEEK_SET)) == -1); k++) usleep (50000); + if (status == -1) { + /* fprintf (stderr, "ERROR: serious problem seeking position: %d\n", errno); */ + /* clear lock at this point? */ + return (0); + } + } + return (1); +} Index: /trunk/Ohana/src/skycalc/src/moon.c =================================================================== --- /trunk/Ohana/src/skycalc/src/moon.c (revision 2496) +++ /trunk/Ohana/src/skycalc/src/moon.c (revision 2496) @@ -0,0 +1,465 @@ +# include + +void lpmoon(double jd, double lat, double sid, + double* ra, double* dec, double* dist) { + + +/* implements "low precision" moon algorithms from + Astronomical Almanac (p. D46 in 1992 version). Does + apply the topocentric correction. + Units are as follows + jd,lat, sid; decimal hours + *ra, *dec, decimal hours, degrees + *dist; earth radii */ + + + double T, lambda, beta, pie, l, m, n, x, y, z, alpha, delta, + rad_lat, rad_lst, distance, topo_dist; + char dummy[40]; /* to fix compiler bug on IBM system */ + + T = (jd - J2000) / 36525.; /* jul cent. since J2000.0 */ + + lambda = 218.32 + 481267.883 * T + + 6.29 * sin((134.9 + 477198.85 * T) / DEG_IN_RADIAN) + - 1.27 * sin((259.2 - 413335.38 * T) / DEG_IN_RADIAN) + + 0.66 * sin((235.7 + 890534.23 * T) / DEG_IN_RADIAN) + + 0.21 * sin((269.9 + 954397.70 * T) / DEG_IN_RADIAN) + - 0.19 * sin((357.5 + 35999.05 * T) / DEG_IN_RADIAN) + - 0.11 * sin((186.6 + 966404.05 * T) / DEG_IN_RADIAN); + lambda = lambda / DEG_IN_RADIAN; + beta = 5.13 * sin((93.3 + 483202.03 * T) / DEG_IN_RADIAN) + + 0.28 * sin((228.2 + 960400.87 * T) / DEG_IN_RADIAN) + - 0.28 * sin((318.3 + 6003.18 * T) / DEG_IN_RADIAN) + - 0.17 * sin((217.6 - 407332.20 * T) / DEG_IN_RADIAN); + beta = beta / DEG_IN_RADIAN; + pie = 0.9508 + + 0.0518 * cos((134.9 + 477198.85 * T) / DEG_IN_RADIAN) + + 0.0095 * cos((259.2 - 413335.38 * T) / DEG_IN_RADIAN) + + 0.0078 * cos((235.7 + 890534.23 * T) / DEG_IN_RADIAN) + + 0.0028 * cos((269.9 + 954397.70 * T) / DEG_IN_RADIAN); + pie = pie / DEG_IN_RADIAN; + distance = 1 / sin(pie); + + l = cos(beta) * cos(lambda); + m = 0.9175 * cos(beta) * sin(lambda) - 0.3978 * sin(beta); + n = 0.3978 * cos(beta) * sin(lambda) + 0.9175 * sin(beta); + + x = l * distance; + y = m * distance; + z = n * distance; /* for topocentric correction */ + + + /* lat isn't passed right on some IBM systems unless you do this + or something like it! */ + sprintf(dummy,"%f",lat); + + rad_lat = lat / DEG_IN_RADIAN; + rad_lst = sid / HRS_IN_RADIAN; + + x = x - cos(rad_lat) * cos(rad_lst); + y = y - cos(rad_lat) * sin(rad_lst); + z = z - sin(rad_lat); + + + topo_dist = sqrt(x * x + y * y + z * z); + + l = x / topo_dist; + m = y / topo_dist; + n = z / topo_dist; + + alpha = atan_circ(l,m); + delta = asin(n); + + *ra = alpha * HRS_IN_RADIAN; + + *dec = delta * DEG_IN_RADIAN; + + *dist = topo_dist; + +} + +/* More accurate (but more elaborate and slower) lunar + ephemeris, from Jean Meeus' *Astronomical Formulae For Calculators*, + pub. Willman-Bell. Includes all the terms given there. */ + +void accumoon (jd,geolat,lst,elevsea,geora,geodec,geodist,topora,topodec,topodist) + double jd,geolat,lst,elevsea; + double *geora,*geodec,*geodist,*topora,*topodec,*topodist; +{ + double pie, dist; /* horiz parallax */ + double Lpr,M,Mpr,D,F,Om,T,Tsq,Tcb; + double e,lambda,B,beta,om1,om2; + double sinx, x, y, z, l, m, n; + double x_geo, y_geo, z_geo; /* geocentric position of *observer* */ + + jd = jd + etcorr(jd)/SEC_IN_DAY; /* approximate correction to ephemeris time */ + T = (jd - 2415020.) / 36525.; /* this based around 1900 ... */ + Tsq = T * T; + Tcb = Tsq * T; + + Lpr = 270.434164 + 481267.8831 * T - 0.001133 * Tsq + + 0.0000019 * Tcb; + M = 358.475833 + 35999.0498*T - 0.000150*Tsq + - 0.0000033*Tcb; + Mpr = 296.104608 + 477198.8491*T + 0.009192*Tsq + + 0.0000144*Tcb; + D = 350.737486 + 445267.1142*T - 0.001436 * Tsq + + 0.0000019*Tcb; + F = 11.250889 + 483202.0251*T -0.003211 * Tsq + - 0.0000003*Tcb; + Om = 259.183275 - 1934.1420*T + 0.002078*Tsq + + 0.0000022*Tcb; + + Lpr = circulo(Lpr); + Mpr = circulo(Mpr); + M = circulo(M); + D = circulo(D); + F = circulo(F); + Om = circulo(Om); + + + sinx = sin((51.2 + 20.2 * T)/DEG_IN_RADIAN); + Lpr = Lpr + 0.000233 * sinx; + M = M - 0.001778 * sinx; + Mpr = Mpr + 0.000817 * sinx; + D = D + 0.002011 * sinx; + + sinx = 0.003964 * sin((346.560+132.870*T -0.0091731*Tsq)/DEG_IN_RADIAN); + + Lpr = Lpr + sinx; + Mpr = Mpr + sinx; + D = D + sinx; + F = F + sinx; + + sinx = sin(Om/DEG_IN_RADIAN); + Lpr = Lpr + 0.001964 * sinx; + Mpr = Mpr + 0.002541 * sinx; + D = D + 0.001964 * sinx; + F = F - 0.024691 * sinx; + F = F - 0.004328 * sin((Om + 275.05 -2.30*T)/DEG_IN_RADIAN); + + e = 1 - 0.002495 * T - 0.00000752 * Tsq; + + M = M / DEG_IN_RADIAN; /* these will all be arguments ... */ + Mpr = Mpr / DEG_IN_RADIAN; + D = D / DEG_IN_RADIAN; + F = F / DEG_IN_RADIAN; + + lambda = Lpr + 6.288750 * sin(Mpr) + + 1.274018 * sin(2*D - Mpr) + + 0.658309 * sin(2*D) + + 0.213616 * sin(2*Mpr) + - e * 0.185596 * sin(M) + - 0.114336 * sin(2*F) + + 0.058793 * sin(2*D - 2*Mpr) + + e * 0.057212 * sin(2*D - M - Mpr) + + 0.053320 * sin(2*D + Mpr) + + e * 0.045874 * sin(2*D - M) + + e * 0.041024 * sin(Mpr - M) + - 0.034718 * sin(D) + - e * 0.030465 * sin(M+Mpr) + + 0.015326 * sin(2*D - 2*F) + - 0.012528 * sin(2*F + Mpr) + - 0.010980 * sin(2*F - Mpr) + + 0.010674 * sin(4*D - Mpr) + + 0.010034 * sin(3*Mpr) + + 0.008548 * sin(4*D - 2*Mpr) + - e * 0.007910 * sin(M - Mpr + 2*D) + - e * 0.006783 * sin(2*D + M) + + 0.005162 * sin(Mpr - D); + + /* And furthermore.....*/ + + lambda = lambda + e * 0.005000 * sin(M + D) + + e * 0.004049 * sin(Mpr - M + 2*D) + + 0.003996 * sin(2*Mpr + 2*D) + + 0.003862 * sin(4*D) + + 0.003665 * sin(2*D - 3*Mpr) + + e * 0.002695 * sin(2*Mpr - M) + + 0.002602 * sin(Mpr - 2*F - 2*D) + + e * 0.002396 * sin(2*D - M - 2*Mpr) + - 0.002349 * sin(Mpr + D) + + e * e * 0.002249 * sin(2*D - 2*M) + - e * 0.002125 * sin(2*Mpr + M) + - e * e * 0.002079 * sin(2*M) + + e * e * 0.002059 * sin(2*D - Mpr - 2*M) + - 0.001773 * sin(Mpr + 2*D - 2*F) + - 0.001595 * sin(2*F + 2*D) + + e * 0.001220 * sin(4*D - M - Mpr) + - 0.001110 * sin(2*Mpr + 2*F) + + 0.000892 * sin(Mpr - 3*D) + - e * 0.000811 * sin(M + Mpr + 2*D) + + e * 0.000761 * sin(4*D - M - 2*Mpr) + + e * e * 0.000717 * sin(Mpr - 2*M) + + e * e * 0.000704 * sin(Mpr - 2 * M - 2*D) + + e * 0.000693 * sin(M - 2*Mpr + 2*D) + + e * 0.000598 * sin(2*D - M - 2*F) + + 0.000550 * sin(Mpr + 4*D) + + 0.000538 * sin(4*Mpr) + + e * 0.000521 * sin(4*D - M) + + 0.000486 * sin(2*Mpr - D); + +/* *eclongit = lambda; */ + + B = 5.128189 * sin(F) + + 0.280606 * sin(Mpr + F) + + 0.277693 * sin(Mpr - F) + + 0.173238 * sin(2*D - F) + + 0.055413 * sin(2*D + F - Mpr) + + 0.046272 * sin(2*D - F - Mpr) + + 0.032573 * sin(2*D + F) + + 0.017198 * sin(2*Mpr + F) + + 0.009267 * sin(2*D + Mpr - F) + + 0.008823 * sin(2*Mpr - F) + + e * 0.008247 * sin(2*D - M - F) + + 0.004323 * sin(2*D - F - 2*Mpr) + + 0.004200 * sin(2*D + F + Mpr) + + e * 0.003372 * sin(F - M - 2*D) + + 0.002472 * sin(2*D + F - M - Mpr) + + e * 0.002222 * sin(2*D + F - M) + + e * 0.002072 * sin(2*D - F - M - Mpr) + + e * 0.001877 * sin(F - M + Mpr) + + 0.001828 * sin(4*D - F - Mpr) + - e * 0.001803 * sin(F + M) + - 0.001750 * sin(3*F) + + e * 0.001570 * sin(Mpr - M - F) + - 0.001487 * sin(F + D) + - e * 0.001481 * sin(F + M + Mpr) + + e * 0.001417 * sin(F - M - Mpr) + + e * 0.001350 * sin(F - M) + + 0.001330 * sin(F - D) + + 0.001106 * sin(F + 3*Mpr) + + 0.001020 * sin(4*D - F) + + 0.000833 * sin(F + 4*D - Mpr); + /* not only that, but */ + B = B + 0.000781 * sin(Mpr - 3*F) + + 0.000670 * sin(F + 4*D - 2*Mpr) + + 0.000606 * sin(2*D - 3*F) + + 0.000597 * sin(2*D + 2*Mpr - F) + + e * 0.000492 * sin(2*D + Mpr - M - F) + + 0.000450 * sin(2*Mpr - F - 2*D) + + 0.000439 * sin(3*Mpr - F) + + 0.000423 * sin(F + 2*D + 2*Mpr) + + 0.000422 * sin(2*D - F - 3*Mpr) + - e * 0.000367 * sin(M + F + 2*D - Mpr) + - e * 0.000353 * sin(M + F + 2*D) + + 0.000331 * sin(F + 4*D) + + e * 0.000317 * sin(2*D + F - M + Mpr) + + e * e * 0.000306 * sin(2*D - 2*M - F) + - 0.000283 * sin(Mpr + 3*F); + + om1 = 0.0004664 * cos(Om/DEG_IN_RADIAN); + om2 = 0.0000754 * cos((Om + 275.05 - 2.30*T)/DEG_IN_RADIAN); + + beta = B * (1. - om1 - om2); + /* *eclatit = beta; */ + + pie = 0.950724 + + 0.051818 * cos(Mpr) + + 0.009531 * cos(2*D - Mpr) + + 0.007843 * cos(2*D) + + 0.002824 * cos(2*Mpr) + + 0.000857 * cos(2*D + Mpr) + + e * 0.000533 * cos(2*D - M) + + e * 0.000401 * cos(2*D - M - Mpr) + + e * 0.000320 * cos(Mpr - M) + - 0.000271 * cos(D) + - e * 0.000264 * cos(M + Mpr) + - 0.000198 * cos(2*F - Mpr) + + 0.000173 * cos(3*Mpr) + + 0.000167 * cos(4*D - Mpr) + - e * 0.000111 * cos(M) + + 0.000103 * cos(4*D - 2*Mpr) + - 0.000084 * cos(2*Mpr - 2*D) + - e * 0.000083 * cos(2*D + M) + + 0.000079 * cos(2*D + 2*Mpr) + + 0.000072 * cos(4*D) + + e * 0.000064 * cos(2*D - M + Mpr) + - e * 0.000063 * cos(2*D + M - Mpr) + + e * 0.000041 * cos(M + D) + + e * 0.000035 * cos(2*Mpr - M) + - 0.000033 * cos(3*Mpr - 2*D) + - 0.000030 * cos(Mpr + D) + - 0.000029 * cos(2*F - 2*D) + - e * 0.000029 * cos(2*Mpr + M) + + e * e * 0.000026 * cos(2*D - 2*M) + - 0.000023 * cos(2*F - 2*D + Mpr) + + e * 0.000019 * cos(4*D - M - Mpr); + + beta = beta/DEG_IN_RADIAN; + lambda = lambda/DEG_IN_RADIAN; + l = cos(lambda) * cos(beta); + m = sin(lambda) * cos(beta); + n = sin(beta); + eclrot(jd,&l,&m,&n); + + dist = 1/sin((pie)/DEG_IN_RADIAN); + x = l * dist; + y = m * dist; + z = n * dist; + + *geora = atan_circ(l,m) * HRS_IN_RADIAN; + *geodec = asin(n) * DEG_IN_RADIAN; + *geodist = dist; + + geocent (lst, geolat, elevsea, &x_geo, &y_geo, &z_geo); + + x = x - x_geo; /* topocentric correction using elliptical earth fig. */ + y = y - y_geo; + z = z - z_geo; + + *topodist = sqrt(x*x + y*y + z*z); + + l = x / (*topodist); + m = y / (*topodist); + n = z / (*topodist); + + *topora = atan_circ(l,m) * HRS_IN_RADIAN; + *topodec = asin(n) * DEG_IN_RADIAN; + +} + +/* returns jd at which moon is at a given + altitude, given jdguess as a starting point. In current version + uses high-precision moon -- execution time does not seem to be + excessive on modern hardware. If it's a problem on your machine, + you can replace calls to 'accumoon' with 'lpmoon' and remove + the 'elevsea' argument. */ +double jd_moon_alt (double alt, double jdguess, double lat, double longit, double elevsea) { + + double jdout; + double deriv, err, del = 0.002; + double ra,dec,dist,geora,geodec,geodist,sid,ha,alt2,alt3,az; + short i = 0; + + /* first guess */ + + sid=lst(jdguess,longit); + accumoon(jdguess,lat,sid,elevsea,&geora,&geodec,&geodist, + &ra,&dec,&dist); + ha = lst(jdguess,longit) - ra; + alt2 = altit(dec,ha,lat,&az); + jdguess = jdguess + del; + sid = lst(jdguess,longit); + accumoon(jdguess,lat,sid,elevsea,&geora,&geodec,&geodist, + &ra,&dec,&dist); + alt3 = altit(dec,(sid - ra),lat,&az); + err = alt3 - alt; + deriv = (alt3 - alt2) / del; + while((fabs(err) > 0.1) && (i < 10)) { + jdguess = jdguess - err/deriv; + sid=lst(jdguess,longit); + accumoon(jdguess,lat,sid,elevsea,&geora,&geodec,&geodist, + &ra,&dec,&dist); + alt3 = altit(dec,(sid - ra),lat,&az); + err = alt3 - alt; + i++; + } + if(i >= 9) jdguess = -1000.; + jdout = jdguess; + return(jdout); +} + +/* Given site position, return Moonset for closest midnight */ +double moonset_tonight (struct SC_date_time date, double lat, double longit, double elevsea, double elev) { + + double jd, jdmid, stmid; + double min_alt, max_alt; + double geora, geodec, geodist; /* geocent for moon, not used here.*/ + double ramoon, decmoon, distmoon; + double hamoonrise, hamoonset, tmoonrise, tmoonset, jdmoonrise, jdmoonset; + short dow; /* day of week */ + struct SC_date_time date_midnight, tmpdate; + double dt, lst0, lst1, djd, horiz; + + horiz = sqrt (2. * elev / 6378140.) * DEG_IN_RADIAN; + + /* find offset in hours from longit to greenwich */ + jd = date_to_jd (date); /* true jd now */ + lst0 = lst (jd, 0.0); /* lst at long = 0 */ + lst1 = lst (jd, longit); /* local lst now */ + dt = lst0 - lst1; + if (dt < 0) dt += 24; + + /* midnight at greenwich */ + date_midnight = date; + date_midnight.h = 0; + date_midnight.mn = 0; + date_midnight.s = 0; + + /* find jd for local midnight, select the *closest* midnight */ + jdmid = date_to_jd (date_midnight) - dt / 24.0; + djd = jd - jdmid; + if (djd < -0.5) jdmid -= 1.0; + if (djd > 0.5) jdmid += 1.0; + stmid = lst (jdmid,longit); + + accumoon (jdmid,lat,stmid,elevsea,&geora,&geodec,&geodist,&ramoon,&decmoon,&distmoon); + + min_max_alt (lat,decmoon,&min_alt,&max_alt); /* rough check -- occurs? */ + if (max_alt < -(0.83+horiz)) return (-1); + if (min_alt > -(0.83+horiz)) return (-1); + + /* compute moonrise and set if they're likely to occur */ + + hamoonset = ha_alt(decmoon,lat,-(0.83+horiz)); /* rough approx. */ + + tmoonset = adj_time(ramoon+hamoonset-stmid); + jdmoonset = jdmid + tmoonset / 24.; + jdmoonset = jd_moon_alt(-(0.83+horiz),jdmoonset,lat,longit,elevsea); + + return (jdmoonset); + +} + +/* Given site position, return Moonrise for closest midnight */ +double moonrise_tonight (struct SC_date_time date, double lat, double longit, double elevsea, double elev) { + + double jd, jdmid, stmid; + double min_alt, max_alt; + double geora, geodec, geodist; /* geocent for moon, not used here.*/ + double ramoon, decmoon, distmoon; + double hamoonrise, hamoonset, tmoonrise, tmoonset, jdmoonrise, jdmoonset; + short dow; /* day of week */ + struct SC_date_time date_midnight, tmpdate; + double dt, lst0, lst1, djd, horiz; + + horiz = sqrt (2. * elev / 6378140.) * DEG_IN_RADIAN; + + /* find offset in hours from longit to greenwich */ + jd = date_to_jd (date); /* true jd now */ + lst0 = lst (jd, 0.0); /* lst at long = 0 */ + lst1 = lst (jd, longit); /* local lst now */ + dt = lst0 - lst1; + if (dt < 0) dt += 24; + + /* midnight at greenwich */ + date_midnight = date; + date_midnight.h = 0; + date_midnight.mn = 0; + date_midnight.s = 0; + + /* find jd for local midnight, select the *closest* midnight */ + jdmid = date_to_jd (date_midnight) - dt / 24.0; + djd = jd - jdmid; + if (djd < -0.5) jdmid -= 1.0; + if (djd > 0.5) jdmid += 1.0; + stmid = lst (jdmid,longit); + + accumoon (jdmid,lat,stmid,elevsea,&geora,&geodec,&geodist,&ramoon,&decmoon,&distmoon); + + min_max_alt (lat,decmoon,&min_alt,&max_alt); /* rough check -- occurs? */ + if (max_alt < -(0.83+horiz)) return (-1); + if (min_alt > -(0.83+horiz)) return (-1); + + /* compute moonrise and set if they're likely to occur */ + + hamoonset = ha_alt(decmoon,lat,-(0.83+horiz)); /* rough approx. */ + + tmoonrise = adj_time(ramoon-hamoonset-stmid); + jdmoonrise = jdmid + tmoonrise / 24.; + jdmoonrise = jd_moon_alt(-(0.83+horiz),jdmoonrise,lat,longit,elevsea); + + return (jdmoonrise); + +} Index: /trunk/Ohana/src/skycalc/src/moon_interference.c =================================================================== --- /trunk/Ohana/src/skycalc/src/moon_interference.c (revision 2496) +++ /trunk/Ohana/src/skycalc/src/moon_interference.c (revision 2496) @@ -0,0 +1,875 @@ +/* + This program segment is extracted from the skycalc program written + by John Thorstensen from Dartmouth College. For question please + contact: John.Thorstensen@Dartmouth.edu + This is one long and ugly, but accurate program. Ideally I should + have split this into separate modules, but this may be replaced by + Bernt's algorithm anyway... so "let's wait and see what happens" + -Rosemary Alles 04/18/2001 + */ + +#include +#include +#include +#include +#include +#include + +/* some (not all) physical, mathematical, and astronomical constants + used are defined here. */ +#define DEG_IN_RADIAN 57.2957795130823 +#define HRS_IN_RADIAN 3.819718634205 +#define EQUAT_RAD 6378137. /* equatorial radius of earth, meters */ +#define TWOPI 6.28318530717959 +#define J2000 2451545. /* Julian date at standard epoch */ +#define SEC_IN_DAY 86400. + +struct coord { + short sign; /* carry sign explicitly since -0 not neg. */ + double hh; + double mm; + double ss; +}; + +struct date_time { + short y; + short mo; + short d; + short h; + short mn; + float s; +}; + +/* Gloabals */ +int update_on = 0; +double update_delta = 0.; + +/* Prototypes */ +void load_site(double* longit, double* lat, double* stdz, + short* use_dst, char* zone_name, char* zabr, + double* elevsea, double* elev, double* horiz, + char* site_name); +double zone(short use_dst, double stdz, + double jd, double jdb, double jde); +void caldat(double jdin, struct date_time* date, short *dow); +short day_of_week(double jd); +double day_of_year(double jd); +void find_dst_bounds(short yr, double stdz, short use_dst, + double *jdb, double* jde); +int get_sys_date(struct date_time *date, short use_dst, short enter_ut, + short night_date, double stdz, double toffset, double* jd); +double lst(double jd, double longit); +double date_to_jd(struct date_time date); +void lpmoon(double jd, double lat, double sid, + double* ra, double* dec, double* dist); +void put_coords(double deci, int precision, int showsign); +double get_coord(void); +void dec_to_bab (double deci, struct coord* bab); +double atan_circ(double x, double y); +double adj_time(double x); + + +void +load_site(double* longit, double* lat, double* stdz, + short* use_dst, char* zone_name, char* zabr, + double* elevsea, double* elev, double* horiz, + char* site_name) { + +/* sets the site-specific quantities; these are + longit = W longitude in decimal hours + lat = N latitude in decimal degrees + stdz = standard time zone offset, hours + elevsea = elevation above sea level (for absolute location) + elev = observatory elevation above horizon, meters + horiz = (derived) added zenith distance for rise/set due + to elevation + use_dst = 0 don't use it + 1 use USA convention + 2 use Spanish convention + < 0 Southern hemisphere (reserved, unimplimented) + zone_name = name of time zone, e. g. Eastern + zabr = single-character abbreviation of time zone + site_name = name of site. */ + + strcpy(site_name, "Mauna Kea, Hawaii"); + strcpy(zone_name, "Hawaiian"); + *zabr = 'H'; + *use_dst = 0; + *longit = 10.36478; + *lat = 19.8267; + *elevsea = 4215.; + *elev = 4215.; /* yow! */ + *stdz = 10.; + + /* now compute derived quantity "horiz" = depression of horizon.*/ + *horiz = sqrt(2. * *elev / EQUAT_RAD) * DEG_IN_RADIAN; +} + +double +zone(short use_dst, double stdz, + double jd, double jdb, double jde) { + + /* Returns zone time offset when standard time zone is stdz, + when daylight time begins (for the year) on jdb, and ends + (for the year) on jde. This is parochial to the northern + hemisphere. */ + /* Extension -- specifying a negative value of use_dst reverses + the logic for the Southern hemisphere; then DST is assumed for + the Southern hemisphere summer (which is the end and beginning + of the year. */ + + if(use_dst == 0) return(stdz); + else if((jd > jdb) && (jd < jde) && (use_dst > 0)) return(stdz-1.); + /* next line .. use_dst < 0 .. for Southern Hemisphere sites. */ + else if(((jd < jdb) || (jd > jde)) && (use_dst < 0)) return(stdz-1.); + else return(stdz); +} + +void +caldat(double jdin, struct date_time* date, short *dow) { + + /* from Jean Meeus, Astronomical Formulae for Calculators, + published by Willman-Bell Inc. + Avoids a copyrighted routine from Numerical Recipes. + Tested and works properly from the beginning of the + Gregorian calendar era (1583) to beyond 3000 AD. */ + + double jdtmp; + long alpha; + long Z; + long A, B, C, D, E; + double F; + double x; /* for day-of-week calculation */ + + jdtmp = jdin + 0.5; + + Z = (long) jdtmp; + + x = Z/7.+0.01; + *dow = 7.*(x - (long) x); /* truncate for day of week */ + + F = jdtmp - Z; + + if(Z < 2299161) A = Z; + else { + alpha = (long) ((Z - 1867216.25) / 36524.25); + A = Z + 1 + alpha - (long) (alpha / 4); + } + + B = A + 1524; + C = ((B - 122.1) / 365.25); + D = (365.25 * C);\ + E = ((B - D) / 30.6001); + + date->d = B - D - (long)(30.6001 * E); + if(E < 13.5) date->mo = E - 1; + else date->mo = E - 13; + if(date->mo > 2.5) date->y = C - 4716; + else date->y = C - 4715; + + date->h = F * 24.; /* truncate */ + date->mn = (F - ((float) date->h)/24.) * 1440.; + date->s = (F - ((float) date->h)/24. - + ((float) date->mn)/1440.) * 86400; + +} + +short +day_of_week(double jd) { + + /* returns day of week for a jd, 0 = Mon, 6 = Sun. */ + + double x; + long i; + short d; + + jd = jd+0.5; + i = jd; /* truncate */ + x = i/7.+0.01; + d = 7.*(x - (long) x); /* truncate */ + return(d); +} + +double +day_of_year(double jd) { + + double jdjan0; + struct date_time date; + short dow; + + caldat(jd,&date,&dow); + /* find jd of "jan 0" = Dec 31 of previous year */ + date.y = date.y - 1; + date.mo = 12; + date.d = 31; + date.h = 0; + date.mn = 0; + date.s = 0.; + jdjan0 = date_to_jd(date); + return(jd - jdjan0); +} + +void +find_dst_bounds(short yr, double stdz, short use_dst, + double *jdb, double* jde) { + + /* finds jd's at which daylight savings time begins + and ends. The parameter use_dst allows for a number + of conventions, namely: + 0 = don't use it at all (standard time all the time) + 1 = use USA convention (1st Sun in April to + last Sun in Oct after 1986; last Sun in April before) + 2 = use Spanish convention (for Canary Islands) + -1 = use Chilean convention (CTIO). + -2 = Australian convention (for AAT). + Negative numbers denote sites in the southern hemisphere, + where jdb and jde are beginning and end of STANDARD time for + the year. + It's assumed that the time changes at 2AM local time; so + when clock is set ahead, time jumps suddenly from 2 to 3, + and when time is set back, the hour from 1 to 2 AM local + time is repeated. This could be changed in code if need be. */ + + struct date_time trial; + + if((use_dst == 1) || (use_dst == 0)) { + /* USA Convention, and including no DST to be defensive */ + /* Note that this ignores various wrinkles such as the + brief Nixon administration flirtation with year-round DST, + the extended DST of WW II, and so on. */ + trial.y = yr; + trial.mo = 4; + if(yr >= 1986) trial.d = 1; + else trial.d = 30; + trial.h = 2; + trial.mn = 0; + trial.s = 0; + + /* Find first Sunday in April for 1986 on ... */ + if(yr >= 1986) + while(day_of_week(date_to_jd(trial)) != 6) + trial.d++; + + /* Find last Sunday in April for pre-1986 .... */ + else while(day_of_week(date_to_jd(trial)) != 6) + trial.d--; + + *jdb = date_to_jd(trial) + stdz/24.; + + /* Find last Sunday in October ... */ + trial.mo = 10; + trial.d = 31; + while(day_of_week(date_to_jd(trial)) != 6) { + trial.d--; + } + *jde = date_to_jd(trial) + (stdz - 1.)/24.; + } + else if (use_dst == 2) { /* Spanish, for Canaries */ + trial.y = yr; + trial.mo = 3; + trial.d = 31; + trial.h = 2; + trial.mn = 0; + trial.s = 0; + + while(day_of_week(date_to_jd(trial)) != 6) { + trial.d--; + } + *jdb = date_to_jd(trial) + stdz/24.; + trial.mo = 9; + trial.d = 30; + while(day_of_week(date_to_jd(trial)) != 6) { + trial.d--; + } + *jde = date_to_jd(trial) + (stdz - 1.)/24.; + } + else if (use_dst == -1) { /* Chilean, for CTIO, etc. */ + /* off daylight 2nd Sun in March, onto daylight 2nd Sun in October */ + trial.y = yr; + trial.mo = 3; + trial.d = 8; /* earliest possible 2nd Sunday */ + trial.h = 2; + trial.mn = 0; + trial.s = 0; + + while(day_of_week(date_to_jd(trial)) != 6) { + trial.d++; + } + *jdb = date_to_jd(trial) + (stdz - 1.)/24.; + /* note jdb is beginning of STANDARD time in south, + hence use stdz - 1. */ + trial.mo = 10; + trial.d = 8; + while(day_of_week(date_to_jd(trial)) != 6) { + trial.d++; + } + *jde = date_to_jd(trial) + stdz /24.; + } + else if (use_dst == -2) { /* For Anglo-Australian Telescope */ + /* off daylight 1st Sun in March, onto daylight last Sun in October */ + trial.y = yr; + trial.mo = 3; + trial.d = 1; /* earliest possible 1st Sunday */ + trial.h = 2; + trial.mn = 0; + trial.s = 0; + + while(day_of_week(date_to_jd(trial)) != 6) { + trial.d++; + } + *jdb = date_to_jd(trial) + (stdz - 1.)/24.; + /* note jdb is beginning of STANDARD time in south, + hence use stdz - 1. */ + trial.mo = 10; + trial.d = 31; + while(day_of_week(date_to_jd(trial)) != 6) { + trial.d--; + } + *jde = date_to_jd(trial) + stdz /24.; + } +} + +int +get_sys_date(struct date_time *date, short use_dst, short enter_ut, + short night_date, double stdz, double toffset, double* jd) { + + /* Reads the system clock; loads up the date structure + to conform to the prevailing conventions for the interpretation + of times. Optionally adds "toffset" minutes to the system + clock, as in x minutes in the future. */ + + time_t t, *tp; + struct tm *stm; + double jdb, jde; + short dow; + + tp = &t; /* have to initialize pointer variable for it to + serve as an argument. */ + + t = time(tp); + if(t == -1) { + printf("error: system time unavailable during calculation of moon interference\n"); + return(-1); + } + stm = localtime(&t); + date->y = (short) (stm->tm_year + 1900); + date->mo = (short) (stm->tm_mon + 1); + date->d = (short) (stm->tm_mday); + date->h = (short) (stm->tm_hour); + date->mn = (short) (stm->tm_min); + date->s = (float) (stm->tm_sec); + + if(toffset != 0.) { + *jd = date_to_jd(*date); + *jd = *jd + toffset / 1440.; + caldat(*jd,date,&dow); + } + + if(enter_ut == 1) { /* adjust if needed */ + find_dst_bounds(date->y,stdz,use_dst,&jdb,&jde); + *jd = date_to_jd(*date); + *jd = *jd + zone(use_dst,stdz,*jd,jdb,jde)/24.; + caldat(*jd,date,&dow); + } + if((night_date == 1) && (date->h < 12)) { + date->d = date->d - 1; + } + + return(0); /* success */ +} + +double +lst(double jd, double longit) { + + /* returns the local MEAN sidereal time (dec hrs) at julian date jd + at west longitude long (decimal hours). Follows + definitions in 1992 Astronomical Almanac, pp. B7 and L2. + Expression for GMST at 0h ut referenced to Aoki et al, A&A 105, + p.359, 1982. On workstations, accuracy (numerical only!) + is about a millisecond in the 1990s. */ + + double t, ut, jdmid, jdint, jdfrac, sid_g; + long jdin, sid_int; + + jdin = jd; /* fossil code from earlier package which + split jd into integer and fractional parts ... */ + jdint = jdin; + jdfrac = jd - jdint; + if(jdfrac < 0.5) { + jdmid = jdint - 0.5; + ut = jdfrac + 0.5; + } + else { + jdmid = jdint + 0.5; + ut = jdfrac - 0.5; + } + t = (jdmid - J2000)/36525; + sid_g = (24110.54841+8640184.812866*t+0.093104*t*t-6.2e-6*t*t*t)/SEC_IN_DAY; + sid_int = sid_g; + sid_g = sid_g - (double) sid_int; + sid_g = sid_g + 1.0027379093 * ut - longit/24.; + sid_int = sid_g; + sid_g = (sid_g - (double) sid_int) * 24.; + if(sid_g < 0.) sid_g = sid_g + 24.; + return(sid_g); +} + + +double +date_to_jd(struct date_time date) { + +/* From Meeus' Astronomical Formulae for Calculators. The two JD + conversion routines routines were replaced 1998 November 29 to + avoid inclusion of copyrighted "Numerical Recipes" code. A test + of 1 million random JDs between 1585 and 3200 AD gave the same + conversions as the NR routines. */ + + + double jd; + int y, m; + long A, B; + + if(date.mo <= 2) { + y = date.y - 1; + m = date.mo + 12; + } + else { + y = date.y; + m = date.mo; + } + + A = (long) (y / 100.); + B = 2 - A + (long) (A / 4.); + + jd = (long) (365.25 * y) + (long) (30.6001 * (m + 1)) + date.d + + 1720994.5; + + jd += date.h / 24. + date.mn / 1440. + date.s / 86400.; + + if(date.y > 1583) return(jd + B); + else return(jd); + /* Not quite right, since Gregorian calendar first + adopted around Oct 1582. But fine for modern. */ +} + +void +lpmoon(double jd, double lat, double sid, + double* ra, double* dec, double* dist) { + + +/* implements "low precision" moon algorithms from + Astronomical Almanac (p. D46 in 1992 version). Does + apply the topocentric correction. + Units are as follows + jd,lat, sid; decimal hours + *ra, *dec, decimal hours, degrees + *dist; earth radii */ + + + double T, lambda, beta, pie, l, m, n, x, y, z, alpha, delta, + rad_lat, rad_lst, distance, topo_dist; + char dummy[40]; /* to fix compiler bug on IBM system */ + + T = (jd - J2000) / 36525.; /* jul cent. since J2000.0 */ + + lambda = 218.32 + 481267.883 * T + + 6.29 * sin((134.9 + 477198.85 * T) / DEG_IN_RADIAN) + - 1.27 * sin((259.2 - 413335.38 * T) / DEG_IN_RADIAN) + + 0.66 * sin((235.7 + 890534.23 * T) / DEG_IN_RADIAN) + + 0.21 * sin((269.9 + 954397.70 * T) / DEG_IN_RADIAN) + - 0.19 * sin((357.5 + 35999.05 * T) / DEG_IN_RADIAN) + - 0.11 * sin((186.6 + 966404.05 * T) / DEG_IN_RADIAN); + lambda = lambda / DEG_IN_RADIAN; + beta = 5.13 * sin((93.3 + 483202.03 * T) / DEG_IN_RADIAN) + + 0.28 * sin((228.2 + 960400.87 * T) / DEG_IN_RADIAN) + - 0.28 * sin((318.3 + 6003.18 * T) / DEG_IN_RADIAN) + - 0.17 * sin((217.6 - 407332.20 * T) / DEG_IN_RADIAN); + beta = beta / DEG_IN_RADIAN; + pie = 0.9508 + + 0.0518 * cos((134.9 + 477198.85 * T) / DEG_IN_RADIAN) + + 0.0095 * cos((259.2 - 413335.38 * T) / DEG_IN_RADIAN) + + 0.0078 * cos((235.7 + 890534.23 * T) / DEG_IN_RADIAN) + + 0.0028 * cos((269.9 + 954397.70 * T) / DEG_IN_RADIAN); + pie = pie / DEG_IN_RADIAN; + distance = 1 / sin(pie); + + l = cos(beta) * cos(lambda); + m = 0.9175 * cos(beta) * sin(lambda) - 0.3978 * sin(beta); + n = 0.3978 * cos(beta) * sin(lambda) + 0.9175 * sin(beta); + + x = l * distance; + y = m * distance; + z = n * distance; /* for topocentric correction */ + + + /* lat isn't passed right on some IBM systems unless you do this + or something like it! */ + sprintf(dummy,"%f",lat); + + rad_lat = lat / DEG_IN_RADIAN; + rad_lst = sid / HRS_IN_RADIAN; + + x = x - cos(rad_lat) * cos(rad_lst); + y = y - cos(rad_lat) * sin(rad_lst); + z = z - sin(rad_lat); + + + topo_dist = sqrt(x * x + y * y + z * z); + + l = x / topo_dist; + m = y / topo_dist; + n = z / topo_dist; + + alpha = atan_circ(l,m); + delta = asin(n); + + *ra = alpha * HRS_IN_RADIAN; + + *dec = delta * DEG_IN_RADIAN; + + *dist = topo_dist; + +} + + +void +put_coords(double deci, int precision, int showsign) { + + +/* prints out a struct coord in a nice format; precision + is a code for how accurate you want it. The options are: + precision = 0; minutes rounded to the nearest minute + precision = 1; minutes rounded to the nearest tenth. + precision = 2; seconds rounded to the nearest second + precision = 3; seconds given to the tenth + precision = 4; seconds given to the hundredth + The program assumes that the line is ready for the coord + to be printed and does NOT deliver a new line at the end + of the output. */ + + + double minutes; /* for rounding off if necess. */ + struct coord out_coord, coords; + char out_string[20]; /* for checking for nasty 60's */ + + dec_to_bab(deci,&coords); /* internally convert to coords*/ + + if(coords.sign == -1) printf("-"); + else printf(" "); /* to preserve alignment */ + + if(precision == 0) { /* round to nearest minute */ + minutes = coords.mm + coords.ss / 60.; + /* check to be sure minutes aren't 60 */ + sprintf(out_string,"%.0f %02.0f",coords.hh,minutes); + sscanf(out_string,"%lf %lf",&out_coord.hh,&out_coord.mm); + if(fabs(out_coord.mm - 60.) < 1.0e-7) { + out_coord.mm = 0.; + out_coord.hh = out_coord.hh + 1.; + } + printf("%2.0f:%02.0f",out_coord.hh,out_coord.mm); + } + + else if(precision == 1) { /* keep nearest tenth of a minute */ + minutes = coords.mm + coords.ss / 60.; + /* check to be sure minutes are not 60 */ + sprintf(out_string,"%.0f %04.1f",coords.hh,minutes); + sscanf(out_string,"%lf %lf",&out_coord.hh, &out_coord.mm); + if(fabs(out_coord.mm - 60.) < 1.0e-7) { + out_coord.mm = 0.; + out_coord.hh = out_coord.hh + 1.; + } + printf("%2.0f:%04.1f", out_coord.hh, out_coord.mm); + } + else if(precision == 2) { + /* check to be sure seconds are not 60 */ + sprintf(out_string,"%.0f %02.0f %02.0f",coords.hh,coords.mm,coords.ss); + sscanf(out_string,"%lf %lf %lf",&out_coord.hh,&out_coord.mm, + &out_coord.ss); + if(fabs(out_coord.ss - 60.) < 1.0e-7) { + out_coord.mm = out_coord.mm + 1.; + out_coord.ss = 0.; + if(fabs(out_coord.mm - 60.) < 1.0e-7) { + out_coord.hh = out_coord.hh + 1.; + out_coord.mm = 0.; + } + } + printf("%2.0f:%02.0f:%02.0f",out_coord.hh,out_coord.mm,out_coord.ss); + } + else if(precision == 3) { + /* the usual shuffle to check for 60's */ + sprintf(out_string,"%.0f %02.0f %04.1f",coords.hh, coords.mm, coords.ss); + sscanf(out_string,"%lf %lf %lf",&out_coord.hh,&out_coord.mm, + &out_coord.ss); + if(fabs(out_coord.ss - 60.) < 1.0e-7) { + out_coord.mm = out_coord.mm + 1.; + out_coord.ss = 0.; + if(fabs(out_coord.mm - 60.) < 1.0e-7) { + out_coord.hh = out_coord.hh + 1.; + out_coord.mm = 0.; + } + } + printf("%2.0f:%02.0f:%04.1f",out_coord.hh,out_coord.mm,out_coord.ss); + } + else { + sprintf(out_string,"%.0f %02.0f %05.2f",coords.hh,coords.mm,coords.ss); + sscanf(out_string,"%lf %lf %lf",&out_coord.hh,&out_coord.mm, + &out_coord.ss); + if(fabs(out_coord.ss - 60.) < 1.0e-6) { + out_coord.mm = out_coord.mm + 1.; + out_coord.ss = 0.; + if(fabs(out_coord.mm - 60.) < 1.0e-6) { + out_coord.hh = out_coord.hh + 1.; + out_coord.mm = 0.; + } + } + printf("%2.0f:%02.0f:%05.2f",out_coord.hh, out_coord.mm, out_coord.ss); + } +} + + +double +get_coord(void) { + +/* Reads a string from the terminal and converts it into + a double-precision coordinate. This is trickier than + it appeared at first, since a -00 tests as non-negative; + the sign has to be picked out and handled explicitly. */ +/* Prompt for input in the calling routine.*/ + + short sign; + double hrs, mins, secs; + char hh_string[6]; /* string with the first coord (hh) */ + char hh1[1]; + short i = 0; + + /* read and handle the hour (or degree) part with sign */ + + scanf("%s",hh_string); + hh1[0] = hh_string[i]; + + while(hh1[0] == ' ') { + /* discard leading blanks */ + i++; + hh1[0] = hh_string[i]; + } + + if(hh1[0] == '-') sign = -1; + + else sign = 1; + + sscanf(hh_string,"%lf", &hrs); + if(sign == -1) hrs = -1. * hrs; + + /* read in the minutes and seconds normally */ + scanf("%lf %lf",&mins,&secs); + + return(sign * (hrs + mins / 60. + secs / 3600.)); +} + +void +dec_to_bab (double deci, struct coord* bab) { + + /* function for converting decimal to babylonian hh mm ss.ss */ + int hr_int, min_int; + + if (deci >= 0.) bab->sign = 1; + else { + bab->sign = -1; + deci = -1. * deci; + } + hr_int = deci; /* use conversion conventions to truncate */ + bab->hh = hr_int; + min_int = 60. * (deci - bab->hh); + bab->mm = min_int; + bab->ss = 3600. * (deci - bab->hh - bab->mm / 60.); +} + +double +atan_circ(double x, double y) { + + /* returns radian angle 0 to 2pi for coords x, y -- + get that quadrant right !! */ + + double theta; + + if((x == 0.) && (y == 0.)) return(0.); /* guard ... */ + + theta = atan2(y,x); /* turns out there is such a thing in math.h */ + while(theta < 0.) theta += TWOPI; + return(theta); +} + + +double +adj_time(double x) { + + /* adjusts a time (decimal hours) to be between -12 and 12, + generally used for hour angles. */ + + if(fabs(x) < 100000.) { /* too inefficient for this! */ + while(x > 12.) { + x = x - 24.; + } + while(x < -12.) { + x = x + 24.; + } + } + + else printf("warning: Out of bounds in adj_time in moon interference routine!\n"); + + return(x); +} + +double +altit(double dec, double ha, double lat, + double* az, double *parang) { + +/* + returns altitude(degr) for dec, ha, lat (decimal degr, hr, degr); + also computes and returns azimuth through pointer argument, + and as an extra added bonus returns parallactic angle (decimal degr) + through another pointer argument. + */ + + double x,y,z; + double sinp, cosp; /* sin and cos of parallactic angle */ + double cosdec, sindec, cosha, sinha, coslat, sinlat; + /* time-savers ... */ + + dec = dec / DEG_IN_RADIAN; + ha = ha / HRS_IN_RADIAN; + lat = lat / DEG_IN_RADIAN; /* thank heavens for pass-by-value */ + cosdec = cos(dec); sindec = sin(dec); + cosha = cos(ha); sinha = sin(ha); + coslat = cos(lat); sinlat = sin(lat); + x = DEG_IN_RADIAN * asin(cosdec*cosha*coslat + sindec*sinlat); + y = sindec*coslat - cosdec*cosha*sinlat; /* due N comp. */ + z = -1. * cosdec*sinha; /* due east comp. */ + *az = atan2(z,y); + + /* as it turns out, having knowledge of the altitude and + azimuth makes the spherical trig of the parallactic angle + less ambiguous ... so do it here! Method uses the + "astronomical triangle" connecting celestial pole, object, + and zenith ... now know all the other sides and angles, + so we can crush it ... */ + + if(cosdec != 0.) { /* protect divide by zero ... */ + sinp = -1. * sin(*az) * coslat / cosdec; + /* spherical law of sines .. note cosdec = sin of codec, + coslat = sin of colat .... */ + cosp = -1. * cos(*az) * cosha - sin(*az) * sinha * sinlat; + /* spherical law of cosines ... also transformed to local + available variables. */ + *parang = atan2(sinp,cosp) * DEG_IN_RADIAN; + /* let the library function find the quadrant ... */ + } + else { /* you're on the pole */ + if(lat >= 0.) *parang = 180.; + else *parang = 0.; + } + + *az *= DEG_IN_RADIAN; /* done with taking trig functions of it ... */ + while(*az < 0.) *az += 360.; /* force 0 -> 360 */ + while(*az >= 360.) *az -= 360.; + + return(x); +} + + +int +main(void) { + + /* Site specific parameters */ + double longit, lat; + double stdz; + short use_dst; + char zabr; + double elevsea; + double elev, horiz; + + char site_name[45]; /* initialized later with + strcpy for portability */ + char zone_name[25]; /* this too */ + + struct date_time date; + double jd, sid; + + /* Position of moon */ + double ra, dec, dist; + double ha, az, par, alt; + + /* Misc stuff */ + short enter_ut = 1; + short night_date = 0; + + /* Load site specific information */ + load_site(&longit,&lat,&stdz,&use_dst,zone_name,&zabr, + &elevsea,&elev,&horiz,site_name); + + /* Get system date and calculate julian date */ + if (get_sys_date(&date, use_dst, enter_ut, night_date, stdz, + update_delta, &jd) != 0) { + printf("error: Can't get system date! \n"); + return(-1); + } + + /* Calcualte local sidereal time */ + sid = lst(jd, longit); + + /* Calculate position of moon */ + lpmoon(jd, lat, sid, &ra, &dec, &dist); + + ha = adj_time(sid - ra); + + alt=altit(dec, ha, lat, &az, &par); + + /* Debug section + printf("moon ra = "); + put_coords(ra, 2, 1); + printf("\n"); + printf("moon dec = "); + put_coords(dec, 2, 1); + printf("\n"); + printf("moon ha = "); + put_coords(ha, 2, 1); + printf("\n"); + printf("moon altitude = %.2f\n", alt); + printf("moon zenith-distance = %.2f\n", 90 - alt); + */ + + /* Calculation of where to move */ + /* If zenith distance is within 40 degrees then we move. + In order to do so the following logic will be applied: + + 1) Is zenith distance within 40 degrees (inclusive) + 2) If so: + if dec > 20 then + move dec by -40 + else + move dec by +40 + + */ + + if (fabs(90.0 - alt) <= 40) { + double adj_dec; + /* Zenith distance within 40 degress */ + if (dec > lat) { + adj_dec = dec - 40.0; + } + else { + adj_dec = dec + 40.0; + } + printf("TRUE "); + put_coords(adj_dec, 2, 1); + printf("\n"); + } + else { + printf("FALSE\n"); + } + + return(0); + +} + + Index: /trunk/Ohana/src/skycalc/src/moondata.c =================================================================== --- /trunk/Ohana/src/skycalc/src/moondata.c (revision 2496) +++ /trunk/Ohana/src/skycalc/src/moondata.c (revision 2496) @@ -0,0 +1,81 @@ +# include +# include +# include + +# define DEG_RAD 57.295779513082322 +# define RAD_DEG 0.017453292519943 + +double atof(); + +void set_site (double *longit, double *lat, double *elevsea, double *elev) { + + *longit = 10.36478; /* W longitude in decimal hours */ + *lat = 19.8267; /* N latitude in decimal degrees */ + *elevsea = 4215.; /* elevation above sea level (for absolute location) */ + *elev = 4215.; /* observatory elevation above horizon, meters */ + +} + +main (int argc, char **argv) { + + struct SC_date_time date, tmpdate; + double longit, lat, elevsea, elev; + double jdnow, ra, dec, sid; + double geora, geodec, geodist, dist; + double RAo, DECo, abx, aby, abz, cs, theta; + double Rsun, Dsun, days; + unsigned long int tzero; + struct tm *stm; + + if (argc != 4) { + fprintf (stderr, "USAGE: moondata (date) (ra) (dec) [ra & dec in dec. deg.]\n"); + exit (1); + } + + RAo = atof (argv[2]); + DECo = atof (argv[3]); + + if (!str_to_time (argv[1], &tzero)) { + fprintf (stderr, "syntax error\n"); + exit (1); + } + stm = gmtime (&tzero); + date.y = (short) (stm->tm_year + 1900); + date.mo = (short) (stm->tm_mon + 1); + date.d = (short) (stm->tm_mday); + date.h = (short) (stm->tm_hour); + date.mn = (short) (stm->tm_min); + date.s = (float) (stm->tm_sec); + + /* + fprintf (stderr, "%4d/%02d/%02d %02d:%02d:%02f\n", date.y, date.mo, date.d, date.h, date.mn, date.s); + */ + + set_site (&longit, &lat, &elevsea, &elev); + jdnow = date_to_jd (date); + /* Calcualte local sidereal time */ + sid = lst(jdnow, longit); + lpsun (jdnow, &Rsun, &Dsun); + + accumoon(jdnow,lat,sid,elevsea,&geora,&geodec,&geodist,&ra,&dec,&dist); + + /* fprintf (stdout, "moon @ %f %f\n", 15*ra, dec); */ + + abx = cos(RAo*RAD_DEG)*cos(DECo*RAD_DEG)*cos(15*ra*RAD_DEG)*cos(dec*RAD_DEG); + aby = sin(RAo*RAD_DEG)*cos(DECo*RAD_DEG)*sin(15*ra*RAD_DEG)*cos(dec*RAD_DEG); + abz = sin(DECo*RAD_DEG)*sin(dec*RAD_DEG); + + cs = abx + aby + abz; + + theta = DEG_RAD * acos (cs); + + days = (Rsun - ra - 12)/24.0; + while (days < -0.5) days += 1.0; + while (days > 0.5) days -= 1.0; + days *= 29.5; + + fprintf (stdout, "moon @ %f %f dist %f %f days from full\n", 15*ra, dec, theta, days); + exit (0); + +} + Index: /trunk/Ohana/src/skycalc/src/skylib.c =================================================================== --- /trunk/Ohana/src/skycalc/src/skylib.c (revision 2496) +++ /trunk/Ohana/src/skycalc/src/skylib.c (revision 2496) @@ -0,0 +1,48 @@ +# include +# include +# include + +void set_site (double *longit, double *lat, double *elevsea, double *elev) { + + *longit = 10.36478; /* W longitude in decimal hours */ + *lat = 19.8267; /* N latitude in decimal degrees */ + *elevsea = 4215.; /* elevation above sea level (for absolute location) */ + *elev = 4215.; /* observatory elevation above horizon, meters */ + +} + +main() + +{ + + struct SC_date_time date, tmpdate; + double longit, lat, elevsea, elev; + double jd; + + set_site (&longit, &lat, &elevsea, &elev); + get_sys_date (&date); + + jd = sunset_tonight (date, lat, longit, elev); + fprintf (stderr, "Sunset (%5.0f m horizon): %f\n", elev, jd); + jd_to_date (jd, &tmpdate); + fprintf (stderr, "%4d/%02d/%02d %02d:%02d:%02f\n", tmpdate.y, tmpdate.mo, tmpdate.d, tmpdate.h, tmpdate.mn, tmpdate.s); + + jd = sunrise_tonight (date, lat, longit, elev); + fprintf (stderr, "Sunrise (%5.0f m horizon): %f\n", elev, jd); + jd_to_date (jd, &tmpdate); + fprintf (stderr, "%4d/%02d/%02d %02d:%02d:%02f\n", tmpdate.y, tmpdate.mo, tmpdate.d, tmpdate.h, tmpdate.mn, tmpdate.s); + + + jd = moonset_tonight (date, lat, longit, elevsea, elev); + fprintf (stderr, "Moonset (%5.0f m horizon): %f\n", elev, jd); + jd_to_date (jd, &tmpdate); + fprintf (stderr, "%4d/%02d/%02d %02d:%02d:%02f\n", tmpdate.y, tmpdate.mo, tmpdate.d, tmpdate.h, tmpdate.mn, tmpdate.s); + + jd = moonrise_tonight (date, lat, longit, elevsea, elev); + fprintf (stderr, "Moonrise (%5.0f m horizon): %f\n", elev, jd); + jd_to_date (jd, &tmpdate); + fprintf (stderr, "%4d/%02d/%02d %02d:%02d:%02f\n", tmpdate.y, tmpdate.mo, tmpdate.d, tmpdate.h, tmpdate.mn, tmpdate.s); + +} + + /* set_zenith (date, lat, longit, objepoch, &objra, &objdec); */ Index: /trunk/Ohana/src/skycalc/src/sun.c =================================================================== --- /trunk/Ohana/src/skycalc/src/sun.c (revision 2496) +++ /trunk/Ohana/src/skycalc/src/sun.c (revision 2496) @@ -0,0 +1,157 @@ +# include + +/* Low precision formulae for the sun, from Almanac p. C24 (1990) */ +/* ra and dec are returned as decimal hours and decimal degrees. */ +void lpsun (double jd, double *ra, double *dec) { + + double n, L, g, lambda,epsilon,alpha,delta,x,y,z; + + n = jd - J2000; + L = 280.460 + 0.9856474 * n; + g = (357.528 + 0.9856003 * n)/DEG_IN_RADIAN; + lambda = (L + 1.915 * sin(g) + 0.020 * sin(2. * g))/DEG_IN_RADIAN; + epsilon = (23.439 - 0.0000004 * n)/DEG_IN_RADIAN; + + x = cos(lambda); + y = cos(epsilon) * sin(lambda); + z = sin(epsilon)*sin(lambda); + + *ra = (atan_circ(x,y))*HRS_IN_RADIAN; + *dec = (asin(z))*DEG_IN_RADIAN; +} + +/* returns jd at which sun is at a given + altitude, given jdguess as a starting point. Uses + low-precision sun, which is plenty good enough. */ +double jd_sun_alt (double alt, double jdguess, double lat, double longit) { + + double jdout; + double deriv, err, del = 0.002; + double ra,dec,ha,alt2,alt3,az; + short i = 0; + + /* first guess */ + + lpsun (jdguess, &ra, &dec); + ha = lst (jdguess,longit) - ra; + alt2 = altit (dec,ha,lat,&az); + jdguess = jdguess + del; + lpsun (jdguess,&ra,&dec); + alt3 = altit(dec,(lst(jdguess,longit) - ra),lat,&az); + err = alt3 - alt; + deriv = (alt3 - alt2) / del; + while((fabs(err) > 0.1) && (i < 10)) { + jdguess = jdguess - err/deriv; + lpsun(jdguess,&ra,&dec); + alt3 = altit(dec,(lst(jdguess,longit) - ra),lat,&az); + err = alt3 - alt; + i++; + if(i == 9) printf ("Sunrise, set, or twilight calculation not converging!\n"); + } + if(i >= 9) jdguess = -1000.; + jdout = jdguess; + return(jdout); +} + +/* Given site position, prints Sun info for the given night. +/* dates are all in UT now */ +double sunset_tonight (struct SC_date_time date, double lat, double longit, double elev) { + + double jd, jdmid0, jdmid, stmid; + double rasun, decsun, min_alt, max_alt, horiz; + double hasunset, jdsunset, jdsunrise, sid; + double dt, lst0, lst1, djd; + struct SC_date_time date_midnight, tmpdate; + short dow; + + horiz = sqrt (2. * elev / 6378140.) * DEG_IN_RADIAN; + + /* find offset in hours from longit to greenwich */ + jd = date_to_jd (date); /* true jd now */ + lst0 = lst (jd, 0.0); /* lst at long = 0 */ + lst1 = lst (jd, longit); /* local lst now */ + dt = lst0 - lst1; + if (dt < 0) dt += 24; + + /* midnight at greenwich */ + date_midnight = date; + date_midnight.h = 0; + date_midnight.mn = 0; + date_midnight.s = 0; + + /* find jd for local midnight, select the *closest* midnight */ + jdmid0 = date_to_jd (date_midnight); + jdmid = jdmid0 + dt / 24.0; + djd = jd - jdmid; + if (djd < -0.5) jdmid -= 1.0; + if (djd > 0.5) jdmid += 1.0; + stmid = lst (jdmid,longit); + + /* sunset / sunrise hour angle */ + lpsun (jdmid, &rasun, &decsun); + hasunset = ha_alt (decsun, lat, -(0.83+horiz)); + if(hasunset > 900.) { /* flag for never sets */ + return (-1); + } + if(hasunset < -900.) { + return (-1); + } + + /* find sunset time */ + jdsunset = jdmid + adj_time(rasun+hasunset-stmid)/24.; + jdsunset = jd_sun_alt (-(0.83+horiz),jdsunset,lat,longit); + + return (jdsunset); + +} + +/* Given site position, prints Sun info for the given night. +/* dates are all in UT now */ +double sunrise_tonight (struct SC_date_time date, double lat, double longit, double elev) { + + double jd, jdmid, stmid; + double rasun, decsun, min_alt, max_alt, horiz; + double hasunset, jdsunset, jdsunrise, sid; + double dt, lst0, lst1, djd; + struct SC_date_time date_midnight, tmpdate; + short dow; + + horiz = sqrt (2. * elev / 6378140.) * DEG_IN_RADIAN; + + /* find offset in hours from longit to greenwich */ + jd = date_to_jd (date); /* true jd now */ + lst0 = lst (jd, 0.0); /* lst at long = 0 */ + lst1 = lst (jd, longit); /* local lst now */ + dt = lst0 - lst1; + if (dt < 0) dt += 24; + + /* midnight at greenwich */ + date_midnight = date; + date_midnight.h = 0; + date_midnight.mn = 0; + date_midnight.s = 0; + + /* find jd for local midnight, select the *closest* midnight */ + jdmid = date_to_jd (date_midnight) - dt / 24.0; + djd = jd - jdmid; + if (djd < -0.5) jdmid -= 1.0; + if (djd > 0.5) jdmid += 1.0; + stmid = lst (jdmid,longit); + + /* sunset / sunrise hour angle */ + lpsun (jdmid, &rasun, &decsun); + hasunset = ha_alt (decsun, lat, -(0.83+horiz)); + if(hasunset > 900.) { /* flag for never sets */ + return (-1); + } + if(hasunset < -900.) { + fprintf (stderr, "Sun down all day!\n"); + return (-1); + } + + /* find sunrise time */ + jdsunrise = jdmid + adj_time(rasun-hasunset-stmid)/24.; + jdsunrise = jd_sun_alt(-(0.83+horiz),jdsunrise,lat,longit); + + return (jdsunrise); +} Index: /trunk/Ohana/src/skycalc/src/sundata.c =================================================================== --- /trunk/Ohana/src/skycalc/src/sundata.c (revision 2496) +++ /trunk/Ohana/src/skycalc/src/sundata.c (revision 2496) @@ -0,0 +1,84 @@ +# include +# include +# include + +# define DEG_RAD 57.295779513082322 +# define RAD_DEG 0.017453292519943 +# define dCOS(A) ((double) cos ((double)RAD_DEG*A)) +# define dSIN(A) ((double) sin ((double)RAD_DEG*A)) + +double atof(); + +void set_site (double *longit, double *lat, double *elevsea, double *elev) { + + *longit = 10.36478; /* W longitude in decimal hours */ + *lat = 19.8267; /* N latitude in decimal degrees */ + *elevsea = 4215.; /* elevation above sea level (for absolute location) */ + *elev = 4215.; /* observatory elevation above horizon, meters */ + +} + +main (int argc, char **argv) { + + struct SC_date_time date, tmpdate; + double longit, lat, elevsea, elev; + double jdnow, ra, dec, sid, alt, az, sind, sinh, cosh; + double geora, geodec, geodist, dist; + double RAo, DECo, abx, aby, abz, cs, theta; + double Rsun, Dsun, days, Hsun; + unsigned long int tzero; + struct tm *stm; + + if (argc != 4) { + fprintf (stderr, "USAGE: sundata (date) (ra) (dec) [ra & dec in dec. deg.]\n"); + exit (1); + } + + RAo = atof (argv[2]); + DECo = atof (argv[3]); + + if (!str_to_time (argv[1], &tzero)) { + fprintf (stderr, "syntax error\n"); + exit (1); + } + stm = gmtime (&tzero); + date.y = (short) (stm->tm_year + 1900); + date.mo = (short) (stm->tm_mon + 1); + date.d = (short) (stm->tm_mday); + date.h = (short) (stm->tm_hour); + date.mn = (short) (stm->tm_min); + date.s = (float) (stm->tm_sec); + + /* + fprintf (stderr, "%4d/%02d/%02d %02d:%02d:%02f\n", date.y, date.mo, date.d, date.h, date.mn, date.s); + */ + + set_site (&longit, &lat, &elevsea, &elev); + jdnow = date_to_jd (date); + + /* Calcualte local sidereal time */ + sid = lst(jdnow, longit); + lpsun (jdnow, &Rsun, &Dsun); + + /* dot product of unit vectors of (RAo,DECo) & (Rsun,Dsun) */ + abx = dCOS(RAo)*dCOS(DECo)*dCOS(15*Rsun)*dCOS(Dsun); + aby = dSIN(RAo)*dCOS(DECo)*dSIN(15*Rsun)*dCOS(Dsun); + abz = dSIN(DECo)*dSIN(Dsun); + cs = abx + aby + abz; + theta = DEG_RAD * acos (cs); + + /***** get sun altitude *****/ + Hsun = 15.0*(sid - Rsun); + + sind = dSIN (Dsun) * dSIN (lat) + dCOS (Dsun) * dCOS (Hsun) * dCOS (lat); + alt = DEG_RAD * asin (sind); + + sinh = - dCOS (Dsun) * dSIN (Hsun); + cosh = dSIN (Dsun) * dCOS (lat) - dCOS (Dsun) * dCOS (Hsun) * dSIN (lat); + az = DEG_RAD * atan2 (sinh, cosh); + + fprintf (stdout, "sun @ %f %f dist %f altaz: %f %f\n", 15*Rsun, Dsun, theta, alt, az); + exit (0); + +} + Index: /trunk/Ohana/src/skycalc/src/time.c =================================================================== --- /trunk/Ohana/src/skycalc/src/time.c (revision 2496) +++ /trunk/Ohana/src/skycalc/src/time.c (revision 2496) @@ -0,0 +1,148 @@ +# include + +/* fill date structure with current date & time (UT) */ +int get_sys_date (struct SC_date_time *date) { + + time_t t; + struct tm *stm; + double jd, jdb, jde; + short dow; + + t = time (0); + if (t == -1) return (-1); + + stm = gmtime (&t); + date->y = (short) (stm->tm_year + 1900); + date->mo = (short) (stm->tm_mon + 1); + date->d = (short) (stm->tm_mday); + date->h = (short) (stm->tm_hour); + date->mn = (short) (stm->tm_min); + date->s = (float) (stm->tm_sec); + + return (0); /* success */ + +} + +/* convert a UT date to JD ( 1900 -- 2100?) */ +double date_to_jd (struct SC_date_time date) { + + short yr1=0, mo1=1; + long jdzpt = 1720982, jdint, inter; + double jd,jdfrac; + + if ((date.y <= 1900) | (date.y >= 2100)) return (0.0); + + if (date.mo <= 2) { + yr1 = -1; + mo1 = 13; + } + + jdint = 365.25*(date.y+yr1); /* truncates */ + inter = 30.6001*(date.mo+mo1); + jdint = jdint+inter+date.d+jdzpt; + jd = jdint; + jdfrac=date.h/24.+date.mn/1440.+date.s/SEC_IN_DAY; + if (jdfrac < 0.5) { + jdint--; + jdfrac=jdfrac+0.5; + } + else jdfrac=jdfrac-0.5; + jd = jdint+jdfrac; + return (jd); +} + +/* convert JD to a UT date & time */ +void jd_to_date (double jdin, struct SC_date_time *date) { + +#define IGREG 2299161 + + /* Adapted from Press, Flannery, Teukolsky, & + Vetterling, Numerical Recipes in C, (Cambridge + University Press), 1st edn, p. 12. */ + + int mm, id, iyyy; /* their notation */ + long ja, jdint, jalpha, jb, jc, jd, je; + float jdfrac; + double x; + + jdin = jdin + 0.5; /* adjust for 1/2 day */ + jdint = jdin; + x = jdint/7.+0.01; + jdfrac = jdin - jdint; + date->h = jdfrac * 24; /* truncate */ + date->mn = (jdfrac - ((float) date->h)/24.) * 1440.; + date->s = (jdfrac - ((float) date->h)/24. - + ((float) date->mn)/1440.) * SEC_IN_DAY; + + if(jdint > IGREG) { + jalpha=((float) (jdint-1867216)-0.25)/36524.25; + ja=jdint+1+jalpha-(long)(0.25*jalpha); + } + else + ja=jdint; + jb=ja+1524; + jc=6680.0+((float) (jb-2439870)-122.1)/365.25; + jd=365*jc+(0.25*jc); + je=(jb-jd)/30.6001; + id=jb-jd-(int) (30.6001*je); + mm=je-1; + if(mm > 12) mm -= 12; + iyyy=jc-4715; + if(mm > 2) --iyyy; + if (iyyy <= 0) --iyyy; + date->y = iyyy; + date->mo = mm; + date->d = id; + +} + + +/* returns the local MEAN sidereal time (dec hrs) at julian date jd + at west longitude long (decimal hours). Follows + definitions in 1992 Astronomical Almanac, pp. B7 and L2. + Expression for GMST at 0h ut referenced to Aoki et al, A&A 105, + p.359, 1982. On workstations, accuracy (numerical only!) + is about a millisecond in the 1990s. */ +double lst (double jd, double longit) { + + double t, ut, jdmid, jdint, jdfrac, sid_g, sid; + long sid_int; + + jdint = (int) jd; + jdfrac = jd - jdint; + + if (jdfrac < 0.5) { + jdmid = jdint - 0.5; + ut = jdfrac + 0.5; + } + else { + jdmid = jdint + 0.5; + ut = jdfrac - 0.5; + } + t = (jdmid - J2000)/36525; + sid_g = (24110.54841+8640184.812866*t+0.093104*t*t-6.2e-6*t*t*t)/SEC_IN_DAY; + sid_int = sid_g; + sid_g = sid_g - (double) sid_int; + sid_g = sid_g + 1.0027379093 * ut - longit/24.; + sid_int = sid_g; + sid_g = (sid_g - (double) sid_int) * 24.; + if (sid_g < 0.) sid_g = sid_g + 24.; + return (sid_g); +} + +/* force time domain to be -12h and 12h */ +double adj_time (double x) { + + /* ridiculously inefficient - use modulo and fractions.. */ + if(fabs(x) < 100000.) { + while(x > 12.) { + x = x - 24.; + } + while(x < -12.) { + x = x + 24.; + } + } + else printf ("Out of bounds in adj_time!\n"); + return(x); +} + Index: /trunk/Ohana/src/skyprobe/skyproberedccd =================================================================== --- /trunk/Ohana/src/skyprobe/skyproberedccd (revision 2496) +++ /trunk/Ohana/src/skyprobe/skyproberedccd (revision 2496) @@ -0,0 +1,44 @@ +#!/bin/sh + +echo " FLIPS - skyproberedccd (imred)" +if [ $# -ne 2 ] ; then + echo " Process: single SkyProbe CCD image reduction with Dark & Flat-Field" + echo " Syntax: skyproberedccd FILE_IN FILE_OUT" + exit 1 +fi + +# Check if INPUT FITS file exists +if [ ! -f $1 ]; then + echo "FITS file '$1' (FILE_IN = 1st argument) does not exist" + exit 1 +fi + +# Create temporary parameter file name +TEMP="/tmp/redccd$$" + +# Prepare parameter file +cat ~skyprobe/.skyprobeimred.param.CCD | sed "s|INFILE|$1|" | sed "s|OUTFILE|$2|" > $TEMP + +# Set a default size for the buffer height +echo "BUFFER" >> $TEMP +echo " height 240" >> $TEMP +echo " end" >> $TEMP +# CPU type: for swap or not +echo "CPU" >> $TEMP +echo " swap 1" >> $TEMP +echo " end" >> $TEMP +# Timing info flag +echo "TIMING" >> $TEMP +echo " timing 1" >> $TEMP +echo " end" >> $TEMP +# End of the parfile +echo "END" >> $TEMP + +# Combine files +imred_lx $TEMP + +# Remove temporary file +rm -f $TEMP +rm -f imred-stat.dat + +exit 0