Author Topic: Sky map and astronomical functions  (Read 12854 times)

Offline z993126

  • BFF
  • ***
  • Posts: 128
    • View Profile
Sky map and astronomical functions
« on: July 23, 2011, 12:54:55 AM »
This is the object that displays a fancy chart of what's in the sky given a particular location (set internally at the moment) for the current time.  Do feel free to laugh or cry at, depending on how the spirit moves you, and feel free to modify it so it makes more sense/is more efficient--just let me see what the modifications were so I can apply them to my own copies.  Or not, it's not like I'll know if you don't.
Code: [Select]
/* ** HEADER for ../obj/sky.c **
2005-Jul-18 - T. Cook
object to display a map of the sky
-----
2011-Jul-21 - T. Cook modified for DeadSouls
2009-Dec-20 - T. Cook modified for Tublib
2005-Oct-04 - T. Cook tidied up
*/
#pragma strong_types

#include <lib.h>
inherit LIB_ITEM;
inherit LIB_ACTIVATE;

#include "../customdefs.h"

#include "../../inc/astro.h"
#include "../../inc/math.h"

inherit MY_FUNS "/math.c";
inherit MY_FUNS "/geo/timezone.c";
inherit MY_FUNS "/astro/chronfuns.c";
inherit MY_FUNS "/astro/earth_sun.c";
inherit MY_FUNS "/astro/earth_moon.c";
inherit MY_FUNS "/astro/planets.c";
inherit MY_FUNS "/astro/stars.c";
inherit MY_FUNS "/geo/regional_divisions.c";
inherit MY_FUNS "/directions.c";
inherit MY_FUNS "/boxdrawing.c";

static void create(){
item::create();
SetKeyName( "skymap" );
SetId( ({ "sky map", "skymap", "map" }) );
SetShort( "a sky map" );
SetLong( "This is an object that will display a map of the sky when activated." );
SetMass( 10 );
SetBaseCost( 0 );
}

void init(){
item::init();
}

int eventTurnOn(){
//-- declare variables
int i_count, i_graphY, i_pos, i_tsx, i_tsy, i_width, i_xf, i_xl, i_yf;
float f_LST, f_xe, f_ye;
string s_abc, s_def;
string *s_ds;
mixed *ma_BodyHA, *ma_BodyAltAz, *ma_BodyRADec, *ma_BP, *ma_JD, *ma_locus, *ma_BodyAbbr, *ma_tD, *ma_tF;
//-----
i_width = this_player()->GetScreen()[0]-3;
if( i_width == 0 ){ i_width = 70; }
switch( lower_case( this_player()->GetName() ) ){
case "ardneh":
ma_locus = ({ -6.0, "US-IA", "Henry County", "Wayland" });
ma_locus += city_loci( ma_locus[1], ma_locus[2], ma_locus[3] );
break;
default:
ma_locus = ({ 0.0, "xx-xx", "", "" });
ma_locus += city_loci( ma_locus[1], ma_locus[2], ma_locus[3] );
break;
}
ma_JD = Utime2JD();
i_pos=0;
// sun, moon, planets, stars
ma_BodyRADec = allocate( 122 );
for( i_count = 0; i_count < 122; ++i_count ){ ma_BodyRADec[i_count] = allocate( 2 ); }
ma_BodyAltAz = allocate( 122 );
for( i_count = 0; i_count < 122; ++i_count ){ ma_BodyAltAz[i_count] = allocate( 2 ); }
ma_BodyHA    = allocate( 122 );
ma_BodyAbbr  = allocate( 122 );
ma_tF        = allocate( 122 );
for( i_count = 0; i_count < 122; ++i_count ){ ma_tF[i_count] = allocate( 2 ); }
f_LST = HHmmss2mm(sprintf("%s",JD2LST(ma_JD[0],ma_JD[1],ma_locus[5])[0..7]));
// luna
ma_BP=calcMoonPos(ma_JD[0],ma_JD[1]);
ma_BodyRADec[0] = ({ ma_BP[0], ma_BP[1] });
ma_BodyHA   [0] = calcHourAngle( f_LST / 60.0, ma_BodyRADec[0][0] );
ma_BodyAltAz[0] = Equatorial2Horizontal(ma_locus[4],ma_BodyRADec[0][1],ma_BodyHA[0]);
switch( MoonPhaseName( ma_JD[0] ) ){
case "full":            ma_BodyAbbr[0] = "%^BOLD%^%^YELLOW%^O%^RESET%^";break;
case "waning gibbous":  ma_BodyAbbr[0] = "%^BOLD%^%^YELLOW%^C%^RESET%^";break;
case "third quarter":   ma_BodyAbbr[0] = "%^BOLD%^%^YELLOW%^C%^RESET%^";break;
case "waning crescent": ma_BodyAbbr[0] = "%^BOLD%^%^YELLOW%^(%^RESET%^";break;
case "new":             ma_BodyAbbr[0] = "%^BOLD%^%^BLACK%^O%^RESET%^";break;
case "waxing crescent": ma_BodyAbbr[0] = "%^BOLD%^%^YELLOW%^)%^RESET%^";break;
case "first quarter":   ma_BodyAbbr[0] = "%^BOLD%^%^YELLOW%^D%^RESET%^";break;
case "waxing gibbous":  ma_BodyAbbr[0] = "%^BOLD%^%^YELLOW%^D%^RESET%^";break;
default:                ma_BodyAbbr[0] = "%^BOLD%^%^YELLOW%^l%^RESET%^";break;
}
// sol
ma_BP=calcSunPos(ma_JD[0],ma_JD[1]);
ma_BodyRADec[1] = ({ ma_BP[0], ma_BP[1] });
ma_BodyHA   [1] = calcHourAngle( f_LST / 60.0, ma_BodyRADec[1][0] );
ma_BodyAltAz[1] = Equatorial2Horizontal(ma_locus[4],ma_BodyRADec[1][1],ma_BodyHA[1]);
ma_BodyAbbr [1] = "%^BOLD%^%^YELLOW%^@%^RESET%^";
// mercury, venus, mars, jupiter, saturn, uranus, neptune,  pluto
for(i_pos=2;i_pos<10;i_pos++){
ma_BP=calcPlanetPos(ma_JD[0],ma_JD[1]);
ma_BodyRADec[i_pos] =({ ma_BP[i_pos-2][0], ma_BP[i_pos-2][1] });
ma_BodyHA   [i_pos] =calcHourAngle(f_LST/60.0,ma_BodyRADec[i_pos][0]);
ma_BodyAltAz[i_pos] =Equatorial2Horizontal(ma_locus[4],ma_BodyRADec[i_pos][1],ma_BodyHA[i_pos]);
switch(i_pos){
case 2: ma_BodyAbbr[i_pos]="%^BOLD%^%^YELLOW%^m%^RESET%^";break;
case 3: ma_BodyAbbr[i_pos]="%^BOLD%^%^YELLOW%^v%^RESET%^";break;
case 4: ma_BodyAbbr[i_pos]="%^BOLD%^%^YELLOW%^a%^RESET%^";break;
case 5: ma_BodyAbbr[i_pos]="%^BOLD%^%^YELLOW%^j%^RESET%^";break;
case 6: ma_BodyAbbr[i_pos]="%^BOLD%^%^YELLOW%^s%^RESET%^";break;
case 7: ma_BodyAbbr[i_pos]="%^BOLD%^%^YELLOW%^u%^RESET%^";break;
case 8: ma_BodyAbbr[i_pos]="%^BOLD%^%^YELLOW%^n%^RESET%^";break;
case 9: ma_BodyAbbr[i_pos]="%^BOLD%^%^YELLOW%^p%^RESET%^";break;
default:ma_BodyAbbr[i_pos]="?";break;
}
}
// assign zodiac positions & convert to horizontal coordinates
ma_BodyRADec[10]=({  3.8546, 20.1499 });ma_BodyAbbr[10]="%^MAGENTA%^A%^RESET%^";  // aries
ma_BodyRADec[11]=({  6.0000, 23.4388 });ma_BodyAbbr[11]="%^MAGENTA%^T%^RESET%^";  // taurus
ma_BodyRADec[12]=({  8.1454, 20.1499 });ma_BodyAbbr[12]="%^MAGENTA%^G%^RESET%^";  // gemini
ma_BodyRADec[13]=({ 10.1393, 11.4718 });ma_BodyAbbr[13]="%^MAGENTA%^C%^RESET%^";  // cancer
ma_BodyRADec[14]=({ 12.0000,  0.0000 });ma_BodyAbbr[14]="%^MAGENTA%^L%^RESET%^";  // leo
ma_BodyRADec[15]=({ 13.8607,-11.4718 });ma_BodyAbbr[15]="%^MAGENTA%^V%^RESET%^";  // virgo
ma_BodyRADec[16]=({ 15.8546,-20.1499 });ma_BodyAbbr[16]="%^MAGENTA%^B%^RESET%^";  // libra
ma_BodyRADec[17]=({ 18.0000,-23.4388 });ma_BodyAbbr[17]="%^MAGENTA%^O%^RESET%^";  // scorpio
ma_BodyRADec[18]=({ 20.1454,-20.1499 });ma_BodyAbbr[18]="%^MAGENTA%^S%^RESET%^";  // sagittarius
ma_BodyRADec[19]=({ 22.1392,-11.4718 });ma_BodyAbbr[19]="%^MAGENTA%^R%^RESET%^";  // capricorn
ma_BodyRADec[20]=({  0.0000,  0.0000 });ma_BodyAbbr[20]="%^MAGENTA%^Q%^RESET%^";  // aquarius
ma_BodyRADec[21]=({  1.8607, 11.4718 });ma_BodyAbbr[21]="%^MAGENTA%^E%^RESET%^";  // pisces
for(i_pos=10;i_pos<22;++i_pos){
ma_BodyHA[i_pos]   =calcHourAngle(f_LST/60.0,ma_BodyRADec[i_pos][0]);
ma_BodyAltAz[i_pos]=Equatorial2Horizontal(ma_locus[4],ma_BodyRADec[i_pos][1],ma_BodyHA[i_pos]);
}
// stars
for(i_pos=22;i_pos<122;i_pos++){
ma_BodyRADec[i_pos]=STARS()[i_pos-22];
ma_BodyHA[i_pos]   =calcHourAngle(f_LST/60.0,ma_BodyRADec[i_pos][0]);
ma_BodyAltAz[i_pos]=Equatorial2Horizontal(ma_locus[4],ma_BodyRADec[i_pos][1],ma_BodyHA[i_pos]);
switch(to_int(STARS()[i_pos-22][2]*100)){
case -150..-101:ma_BodyAbbr[i_pos]="%^BOLD%^%^WHITE%^@%^RESET%^";break;
case -100..-051:ma_BodyAbbr[i_pos]="%^BOLD%^%^WHITE%^O%^RESET%^";break;
case -050..-001:ma_BodyAbbr[i_pos]="%^BOLD%^%^WHITE%^X%^RESET%^";break;
case  000.. 049:ma_BodyAbbr[i_pos]="%^BOLD%^%^WHITE%^*%^RESET%^";break;
case  050.. 099:ma_BodyAbbr[i_pos]="o";break;
case  100.. 149:ma_BodyAbbr[i_pos]="x";break;
case  150.. 199:ma_BodyAbbr[i_pos]="%^BOLD%^%^BLACK%^+%^RESET%^";break;
case  200.. 249:ma_BodyAbbr[i_pos]="%^BOLD%^%^BLACK%^,%^RESET%^";break;
case  250.. 299:ma_BodyAbbr[i_pos]="%^BOLD%^%^BLACK%^'%^RESET%^";break;
case  300.. 349:ma_BodyAbbr[i_pos]="%^BOLD%^%^BLACK%^.%^RESET%^";break;
default:        ma_BodyAbbr[i_pos]="?";break;
}
}
// fill data table with horizontal coordinates
for(i_pos=0;i_pos<sizeof(ma_tF);++i_pos){
ma_tF[i_pos]=({ma_BodyAltAz[i_pos][0],ma_BodyAltAz[i_pos][1]});
}
s_abc="";s_def="";
i_tsy=i_width/4+1;i_tsx=i_width;
if(i_tsy%2==1){i_tsy--;}
ma_tD=allocate(i_tsy);
for(i_count=0;i_count<i_tsy;++i_count){
ma_tD[i_count]=allocate(i_tsx," ");
}
ma_tD[(i_tsy-2)/2]=explode(sprintf("%"+i_tsx+"'-'s",""),"");
// set up direction labels (azimuthal)
s_ds=({"N ","NE","E ","SE","S ","SW","W ","NW"," N"});
i_xl=to_int((i_tsx-2)/8.0);
for(i_pos=0;i_pos<=3;++i_pos){
ma_tD[i_tsy-1][(0+i_xl*i_pos)]=s_ds[i_pos][0..0];
ma_tD[i_tsy-1][(1+i_xl*i_pos)]=s_ds[i_pos][1..1];
ma_tD[i_tsy-1][(i_tsx-2-i_xl*i_pos)]=s_ds[8-i_pos][0..0];
ma_tD[i_tsy-1][(i_tsx-1-i_xl*i_pos)]=s_ds[8-i_pos][1..1];
}
ma_tD[i_tsy-1][(i_tsx/2)]=s_ds[4][0..0];
// fill display with data
i_pos=0;
for(i_pos=sizeof(ma_tF)-1;i_pos>-1;--i_pos){ // backwards to ensure sun/moon visible
f_ye = fwrap( ( to_float( ma_tF[i_pos][0] ) + 90 ) / 180, 0.0, 1.0 ) * ( i_tsy - 3 );
f_xe = fwrap( to_float( ma_tF[i_pos][1] ) / 360, 0.0, 1.0 ) * i_tsx;
if( ma_tF[i_pos][0] != 0 ){
i_yf = i_tsy - to_int( f_ye ) - 3;
i_xf = to_int( f_xe ) + 1;
if( i_yf < 1 ){
write( sprintf( "Error in y value at %d: f_ye=%f  i_yf=%d\n", i_pos, f_ye, i_yf ) );
break;
}
if( i_xf < 1 ){
write( sprintf( "Error in x value at %d: f_xe=%f  i_xf=%d\n", i_pos, f_xe, i_xf ) );
break;
}
ma_tD[i_yf][i_xf-1] = ma_BodyAbbr[i_pos];
}
}
// graph borders
for( i_graphY = 0; i_graphY < sizeof( ma_tD ); ++i_graphY ){
if( i_graphY > 0 && i_graphY < sizeof( ma_tD ) - 2 ){
s_def += "|" + implode( ma_tD[i_graphY], "" ) + "|\n";
}else if( i_graphY == 0 ){
s_def += "+" + implode( explode( sprintf( "%" + i_tsx + "'-'s", "" ), "" ), "" ) + "+\n";
}else if( i_graphY == sizeof( ma_tD ) - 2 ){
s_def += "+" + implode( explode( sprintf( "%" + i_tsx + "'-'s", "" ), "" ), "" ) + "+\n";
}else{
s_def += " " + implode( ma_tD[i_graphY], "" ) + "\n";
}
}
// graph labels
s_abc+=sprintf("%|"+(i_tsx-2)+"s\n","Sky Map for "+UTC2Local(JD2UTC(ma_JD[0],ma_JD[1]),ma_locus[0],ma_locus[1]));
s_abc+=sprintf("%|"+(i_tsx-2)+"s\n",
sprintf(
"at %7.4f degrees %s, %8.4f degrees %s",
abs(ma_locus[4]),
ma_locus[4]>=0.0?"North":"South",
abs(ma_locus[5]),
ma_locus[5]>=0.0?"East":"West"
)
);
s_abc+=s_def;
s_abc+=sprintf(" %|"+i_tsx+"s\n","Direction");
// graph legend
s_abc += "                           Legend:\n";
s_abc += sprintf("%70'-'s\n","");
s_abc += "%^BOLD%^Stars (magnitude)                   Other celestial bodies%^RESET%^\n";
s_abc += "%^BOLD%^%^WHITE%^@%^RESET%^ -1.5 to -1.0    ";
s_abc += "%^BOLD%^%^WHITE%^O%^RESET%^ -1.0 to -0.5    ";
s_abc += ma_BodyAbbr[01] + "  sol            ";
s_abc += ma_BodyAbbr[00] + "  luna (" + MoonPhaseName( ma_JD[0] ) + ")\n";
s_abc += "%^BOLD%^%^WHITE%^X%^RESET%^ -0.5 to  0.0    ";
s_abc += "%^BOLD%^%^WHITE%^*%^RESET%^  0.0 to  0.5    ";
s_abc += ma_BodyAbbr[02] + "  mercury        ";
s_abc += ma_BodyAbbr[03] + "  venus\n";
s_abc += "o  0.5 to  1.0    x  1.0 to  1.5    ";
s_abc += ma_BodyAbbr[04] + "  mars           ";
s_abc += ma_BodyAbbr[05] + "  jupiter\n";
s_abc += "%^BOLD%^%^BLACK%^+%^RESET%^  1.5 to  2.0    ";
s_abc += "%^BOLD%^%^BLACK%^,%^RESET%^  2.0 to  2.5    ";
s_abc += ma_BodyAbbr[06] + "  saturn         ";
s_abc += ma_BodyAbbr[07] + "  uranus\n";
s_abc += "%^BOLD%^%^BLACK%^'%^RESET%^  2.5 to  3.0    ";
s_abc += "%^BOLD%^%^BLACK%^.%^RESET%^  3.0 to  3.5    ";
s_abc += ma_BodyAbbr[08] + "  neptune        ";
s_abc += ma_BodyAbbr[09] + "  pluto\n%^RESET%^";
s_abc += "%^BOLD%^Constellations%^RESET%^\n";
s_abc += ma_BodyAbbr[10] + "  Aries           " + ma_BodyAbbr[16] + "  Libra\n";
s_abc += ma_BodyAbbr[11] + "  Taurus          " + ma_BodyAbbr[17] + "  Scorpio\n";
s_abc += ma_BodyAbbr[12] + "  Gemini          " + ma_BodyAbbr[18] + "  Sagittarius\n";
s_abc += ma_BodyAbbr[13] + "  Cancer          " + ma_BodyAbbr[19] + "  Capricorn\n";
s_abc += ma_BodyAbbr[14] + "  Leo             " + ma_BodyAbbr[20] + "  Aquarius\n";
s_abc += ma_BodyAbbr[15] + "  Virgo           " + ma_BodyAbbr[21] + "  Pisces";
//-----
write( s_abc );
return 1;
}
/* EOF */

What follows are the files it references:
Code: [Select]
/* ** HEADER for ../inc/astro.h **
2009-Mar-15 - T. Cook
Earth-centric astronomical constants
-----
2011-Jul-17 - T. Cook began altering for DeadSouls
2009-Mar-15 - T. Cook separated from math.c
*/
#define Epoch 2444238.5          // 1980 January 0.0
// Constants defining the Sun's apparent orbit, epoch 1980.0
#define Elonge   278.833540      // ecliptic lon. @ epoch
#define Elongp   282.596403      // ecliptic lon. @ perigee
#define Eccent     0.016718      // Earth orbit eccentricity
#define Sunsmax   149598500    // Earth orbit semimaj. axis, km
#define Sunangsiz  0.533128      // sun angular size, deg., @ Sunsmax
// Elements of the Moon's orbit, epoch 1980.0
#define Mmlong    64.975464     // mean lon. of moon @ epoch
#define Mmlongp  349.383063     // mean lon. of perigee @ epoch
#define Mlnode   151.950429     // mean lon. of node @ epoch
#define Minc       5.145396     // moon orbit inclination
#define Mecc       0.054900     // moon orbit eccentricity
#define Mangsiz    0.5181       // angular size @ dist. a from Earth
#define Msmax 384401.0          // moon orbit semimaj. axis, km
#define Mparallax  0.9507       // parallax @ dist. a from Earth
#define Synmonth  29.53058868   // synodic month length (new to new)
/* EOF */

Code: [Select]
/* ** HEADER for ../inc/math.h **
2009-Mar-15 - T. Cook
mathematical constants
-----
2011-Jul-17 - T. Cook began altering for DeadSouls
2009-Mar-15 - T. Cook separated from math.c
*/
#define R12_2 1.05946309435929526456 // 12th root of 2
#define SR2   1.41421356237309504880 // square root of 2
#define PHI   1.61803398874989484820 // golden ratio
#define SR3   1.73205080756887729352 // square root of 3
#define SR5   2.23606797749978969640 // square root of 5
#define DS    2.41421356237309504880 // silver ratio
#define K     2.58498175957925321706 // Sierpinski's constant
#define E     2.71828182845904523536 // natural logarithm
#define PI    3.14159265358979323846 // ratio of circle cicumf. to diam.
/* EOF */

Code: [Select]
/* ** HEADER for ../funs/math.c **
2005-Apr-19 - T. Cook
math functions
-----
2011-Jul-17 - T. Cook began altering for DeadSouls
2009-May-31 - T. Cook added dms2dd(), modified formatting for notepad++
2009-Apr-14 - T. Cook added imod(), fixed iwrap()
2009-Mar-19 - T. Cook tried to fix iwrap()
2009-Mar-15 - T. Cook separated constants to /inc/
2009-Mar-15 - T. Cook added in/cm conversions
2009-Feb-11 - T. Cook modified for tublib
2005-Sep-27 - T. Cook tidied up
to do:
make int-based math funs for semi-arbitrary-precision maths
*/
#pragma strong_types

#include "../inc/math.h"

float rad2deg(mixed angleRad);
float deg2rad(mixed angleDeg);
float fixangle(float angle);
float fwrap(float n,float min,float max);
int iwrap(mixed n,int min,int max);
float *quad(float ym,float yz,float yp);
float calcHourAngle(mixed LST,mixed RA);
float *Equatorial2Horizontal(mixed lat,mixed decl,mixed HA);
float C2F(mixed Celsius);
float F2C(mixed Farenheit);
float cm2in(mixed cm);
float in2cm(mixed inches);
float fmod(mixed a,mixed b);
int imod(int a,int b);
mixed *cut_float(float a);
mixed *dd2dms(float dd);
float dms2dd(int d, int m, mixed s);
float atanb(float y,float x);
/* info for rad2deg()
func: convert radians to degrees
args: mixed angle (radians)
rtrn: float angle (degrees)
*/
float rad2deg(mixed angleRad){
return to_float(180.0 * angleRad / PI);
}
/* info for deg2rad()
func: convert degrees to radians
args: mixed angle (degrees)
rtrn: float angle (radians)
*/
float deg2rad(mixed angleDeg){
return to_float(PI*angleDeg/180.0);
}
/* info for fixangle()
func: convert an arbitrary angle to range 0..359.99
args: float angle
rtrn: float corrected angle
*/
float fixangle(float angle){
return to_float(angle-(360.0*to_int(angle/360.0)));
}
/* info for fwrap()
func: wrap an arbitrary float around specified range
args: float number, lower bound, upper bound
rtrn: float adjusted number
*/
float fwrap(float n,float min,float max){
if(min>=max){return n;}
// this doesn't work with negative n:
return fmod(n-min,(abs(max-min)))+min;
//float range=abs(max-min);
//return min+n-(range*to_int(n/range));  //fixed???
}
/* info for iwrap()
func: wrap an arbitrary int around specified range
args: mixed number, lower bound, upper bound
rtrn: int adjusted number
*/
int iwrap(mixed n,int min,int max){
if(!intp(n)){n=to_int(n);}
if(min>=max){return n;}
return imod((n-min),(abs(max-min)+1))+min;
}
/* info for *quad()
func: finds the parabola throuh the three points (-1,ym),(0,yz),
(1, yp) and returns the coordinates of the max/min (if any)
xe, ye, the values of x where the parabola crosses zero
(roots of the quadratic) and the number of roots (0, 1 or 2)
within the interval [-1, 1]
args: float ym,yz,yp
rtrn: array[nz,z1,z2,x2,ye]
note: "well, this routine is producing sensible answers"
*/
float *quad(float ym,float yz,float yp){
float nz=0,a=0.5*(ym+yp)-yz,b=0.5*(yp-ym),c=yz,xe=-b/(2*a),
ye=(a*xe+b)*xe+c,dis=b*b-4.0*a*c,dx,z1,z2;
if(dis>0){
dx=0.5*sqrt(dis)/abs(a);
z1=xe-dx;
z2=xe+dx;
if(abs(z1)<=1.0) nz+=1;
if(abs(z2)<=1.0) nz+=1;
if(z1<-1.0) z1=z2;
}
return ({nz,z1,z2,xe,ye});
}
/* info for calcHourAngle()
func: calculate hour angle of object
args: float Local Sidereal Time (decimal hours)
      float object Right Ascension (decimal hours)
rtrn: float object hour angle (decimal hours 0..24)
*/
float calcHourAngle(mixed LST,mixed RA){
return fwrap(to_float(LST-RA),0.0,24.0);
}
/* info for *Equatorial2Horizontal()
func: calculate Horizontal coords from Equatorial
args: mixed latitude
      mixed declination of object
      mixed hour angle of object
rtrn: float array ({altitude,azimuth})
*/
float *Equatorial2Horizontal(mixed lat,mixed decl,mixed HA){
float alt, az;
lat=deg2rad(lat);
decl=deg2rad(decl);
HA=deg2rad(HA*15.0);
alt=asin(sin(lat)*sin(decl)+cos(lat)*cos(decl)*cos(HA));
az=acos((cos(lat)*sin(decl)-sin(lat)*cos(decl)*cos(HA))/cos(alt));
az=rad2deg(az);
HA=rad2deg(HA)/15.0;
if(HA>0&&HA<12){az=360-az;}
if(HA>-24&&HA<-12){az=360-az;}
return ({rad2deg(alt),az});
}
/* info for C2F()
func: convert Celsius to Farenheit
args: mixed Celsius
rtrn: float Farenheit
*/
float C2F(mixed Celsius){
return to_float(Celsius*1.8+32);
}
/* info for F2C()
func: convert Farenheit to Celsius
args: mixed Farenheit
rtrn: float Celsius
*/
float F2C(mixed Farenheit){
return to_float((Farenheit-32)/1.8);
}
/* info for cm2in()
func: convert cm to inches
args: mixed cm
rtrn: float inches
*/
float cm2in(mixed cm){
return to_float(cm/2.54);
}
/* info for in2cm()
func: convert inches to cm
args: mixed inches
rtrn: float cm
*/
float in2cm(mixed inches){
return to_float(inches*2.54);
}
/* info for fmod()
func: floating-point-compatible modulus operator
      (sign is same as divisor a la MS Excel)
args: mixed a, b
rtrn: float c
*/
float fmod(mixed a,mixed b){
if(intp(a)){a=to_float(a);}
if(intp(b)){b=to_float(b);}
return abs(a-b*floor((a+0.0)/(b+0.0)))*(abs(b)/b);
}
/* info for imod()
func: integer modulus operator that works like MS Excel's
args: int a, b
rtrn: int c
*/
int imod(int a,int b){
return to_int(abs(a-b*floor((a+0.0)/(b+0.0)))*(abs(b)/b));
}
/* info for *cut_float()
func: chop a floating point number into parts
args: float a
rtrn: mixed array ({int b,decimal d})
*/
mixed *cut_float(float a){
int b;string c;float d;
sscanf(sprintf("%+22.10f",a),"%d.%s",b,c);
d=to_float("0."+c);
return ({b,d});
}
/* info for dd2dms()
func: convert decimal degrees to deg min sec.frac
args: float dd(.frac)
rtrn: mixed array ({dd, mn, sc.frac})
*/
mixed *dd2dms(float dd){
int d=to_int(dd);
int m=abs(to_int((dd-d)*60));
float s=fwrap((((dd-d)*60)-m)*60,0.0,60.0);
//return sprintf("%+04d%s%02d'%07.4f",d,this_player()->query_encoding()==0?"\194\176":" ",m,s);
return ({d,m,s});
}
/* info for dms2dd()
func: convert deg min sec.frac to decimal degrees
args: int deg, int min, mixed sec(.frac)
rtrn: float deg.frac
*/
float dms2dd(int d,int m,mixed s){
return d+m/60.0+s/3600.0;
}
/* info for atanb()
func: atan2
args: float y, float x
rtrn: float atan2
*/
float atanb(float y,float x){
if(y==0.0){return 0.0;}
return 2.0*atan((sqrt(x*x+y*y)-x)/y);
}
/* EOF */

Code: [Select]
/* ** HEADER for ../funs/geo/timezone.c **
2005-Jun-28 - T. Cook
-----
2011-Jul-17 - T. Cook began altering for DeadSouls
2009-Apr-15 - T. Cook modified for tublib
2005-Sep-27 - T. Cook tidied up
*/
#pragma strong_types

inherit "/realms/ardneh/funs/extras.c";

/* info for timezone()
func: get time zone name
args: float time difference from UTC in hours
      string 2-digit country code
      float daylight saving time modifier in hours
rtrn: string array ({abbreviation, zone name })
*/
string *timezone(mixed tz,string cc,float dst){
cc=lower_case(cc[0..1]);
switch(to_int(tz*10)){
case -120: return ({"Y","Yankee Time Zone"});
case -110: return ({"X","X-ray Time Zone"});
case -100: return ({"W","Whiskey Time Zone"});
case -90:
switch(cc){
case "us":case "mx":
if(dst==0.0){ return ({"HAST/AKST","Hawaii-Aleutian/Alaska Standard Time"});}
else{ return ({"HADT/AKDT","Hawaii-Aleutian/Alaska Daylight Time"});}
case "ca":
if(dst==0.0){ return ({"HNY","Heure Normale de Yukon"});}
else{ return ({"HAY","Heure Advance/e de Yukon"});}
        default: return ({"V","Victor Time Zone"});
}
case -80:
switch(cc){
case "us":case "mx":
if(dst==0.0){ return ({"PST","Pacific Standard Time"});}
else{ return ({"PDT","Pacific Daylight Time"});}
case "ca":
if(dst==0.0){ return ({"HNP","Heure Normale du Pacifique"});}
else{ return ({"HAP","Heure Advance/e du Pacifique"});}
default: return ({"U","Uniform Time Zone"});
}
case -70:
switch(cc){
case "us":case "mx":
if(dst==0.0){ return ({"MST","Mountain Standard Time"});}
else{ return ({"MDT","Mountain Daylight Time"});}
case "ca":
if(dst==0.0){ return ({"HNR","Heure Normale des Rocheuses"});}
else{ return ({"HAR","Heure Advance/e des Rocheuses"});}
default: return ({"T","Tango Time Zone"});
}
case -60:
switch(cc){
case "us":case "mx":
if(dst==0.0){ return ({"CST","Central Standard Time"});}
else{ return ({"CDT","Central Daylight Time"});}
case "ca":
if(dst==0.0){ return ({"HNC","Heure Normale du Centre"});}
else{ return ({"HAC","Heure Advance/e du Centre"});}
        default: return ({"S","Sierra Time Zone"});
}
case -50:
switch(cc){
case "us":case "mx":
if(dst==0.0){ return ({"EST","Eastern Standard Time"});}
else{ return ({"EDT","Eastern Daylight Time"});}
case "ca":
if(dst==0.0){ return ({"HNE","Heure Normale de l'Est"});}
else{ return ({"HAE","Heure Advance/e de l'Est"});}
        default: return ({"R","Romeo Time Zone"});
}
case -40:
switch(cc){
case "us":case "mx":
if(dst==0.0){ return ({"AST","Atlantic Standard Time"});}
else{ return ({"ADT","Atlantic Daylight Time"});}
case "ca":
if(dst==0.0){ return ({"HNA","Heure Normale de l'Atlantique"});}
else{ return ({"HAA","Heure Advance/e de l'Atlantique"});}
        default: return ({"Q","Quebec Time Zone"});
}
case -35:
switch(cc){
case "ca":
if(dst==0.0){ return ({"HNT","Heure Normale de Terre-Neuve"});}
else{ return ({"HAT","Heure Advance/e de Terre-Neuve"});}
        default: return ({"",""});
}
case -30: return ({"P","Papa Time Zone"});
    case -20: return ({"O","Oscar Time Zone"});
    case -10: return ({"N","November Time Zone"});
    case 0:
switch(cc){
case "uk":
if(dst==0.0){ return ({"GMT","Greenwich Mean Time"});}
else{ return ({"BST","British Summer Time"});}
case "ie":
if(dst==0.0){ return ({"GMT","Greenwich Mean Time"});}
else{ return ({"IST","Irish Summer Time"});}
case "fo":case "gl":case "is":case "pt":
if(dst==0.0){ return ({"WET","Western European Time"});}
else{ return ({"EST","European Summer Time"});}
        case "bf":case "ci":case "gm":case "gh":case "gw":case "lr":
        case "ml":case "mr":case "ma":case "sh":case "sn":case "sl":
        case "tg": return ({"UTC","Coordinated Universal Time"});
        default: return ({"Z","Zulu Time Zone"});
}
case 10:
switch(cc){
case "dz": return ({"CET","Central European Time"});
case "de":
if(dst==0.0){ return ({"MEZ","Mitteleuropa\"ische Zeit"});}
else{ return ({"MESZ","Mitteleuropa\"ische Sommerzeit"});}
case "al":case "ad":case "at":case "be":case "ba":case "hr":
case "cz":case "dk":case "fx":case "gi":case "hu":case "it":
case "li":case "lu":case "mk":case "mt":case "mc":case "nl":
case "no":case "pl":case "sm":case "sk":case "es":
case "se":case "ch":case "tn":case "va":
if(dst==0.0){ return ({"CET","Central European Time"});}
else{ return ({"CEST","Central European Summer Time"});}
        default: return ({"A","Alpha Time Zone"});
}
case 20:
switch(cc){
case "by":case "bg":case "cy":case "eg":case "ee":case "fi":
case "gr":case "lv":case "md":case "ro":case "ru":case "tr":
case "ua":
if(dst==0.0){ return ({"EET","Eastern European Time"});}
else{ return ({"EEST","Eastern European Summer Time"});}
default: return ({"B","Bravo Time Zone"});
}
case 30:
switch(cc){
case "ru":
if(dst==0.0){ return ({"MSK","Moscow Time"});}
else{ return ({"MSD","Moscow Summer Time"});}
        default: return ({"C","Charlie Time Zone"});
}
case 40: return ({"D","Delta Time Zone"});
case 50: return ({"E","Echo Time Zone"});
case 55:
switch(cc){
case "in": return ({"IST","Indian Standard Time"});
        default: return ({"",""});
}
case 60: return ({"F","Foxtrot Time Zone"});
case 70:
switch(cc){
case "cx": return ({"CXT","Christmas Island Time"});
        default: return ({"G","Golf Time Zone"});
}
case 80:
switch(cc){
case "au": return ({"AWST","Australian Western Standard Time"});
case "cn": return ({"CST","Chinese Standard Time"});
case "hk": return ({"HKT","Hong Kong Time"});
case "sg": return ({"SST","Singapore Standard Time"});
default: return ({"H","Hotel Time Zone"});
}
case 90:
switch(cc){
case "kr": return ({"KST","Korea Standard Time"});
case "jp": return ({"JST","Japan Standard Time"});
default: return ({"I","India Time Zone"});
}
case 95:
switch(cc){
case "au":
if(dst==0.0){ return ({"ACST","Australian Central Standard Time"});}
else{ return ({"ACDT","Australian Central Daylight Time"});}
default: return ({"",""});
}
case 100:
switch(cc){
case "au":
if(dst==0.0){ return ({"ACST","Australian Eastern Standard Time"});}
else{ return ({"AEDT","Australian Eastern Daylight Time"});}
default: return ({"K","Kilo Time Zone"});
}
case 110: return ({"L","Lima Time Zone"});
case 115:
switch(cc){
case "nf": return ({"NFT","Norfolk (Island) Time"});
default: return ({"",""});
}
case 120:
switch(cc){
case "nz":
if(dst==0.0){ return ({"NZST","New Zealand Standard Time"});}
else{ return ({"NZDT","New Zealand Daylight Time"});}
default: return ({"M","Mike Time Zone"});
}
    default: return ({"UTC","Coordinated Universal Time"});
}
}
/* EOF */

Code: [Select]
/* ** HEADER for ../funs/extras.c **
2011-Jul-19 - T. Cook
extra functions that aren't builtin to DeadSouls
*/
#pragma strong_types

string to_string( mixed m_arg );
string indent_string( string s_str );

/* info for to_string()
func: convert a variable to a string
args: mixed argument
rtrn: string argument
*/
string to_string( mixed m_arg ){
if( intp( m_arg ) ){ return sprintf( "%d", m_arg ); }
else if( floatp( m_arg ) ){ return sprintf( "%f", m_arg ); }
else if( stringp( m_arg ) ){ return sprintf( "%s", m_arg ); }
else{ return sprintf( "%O", m_arg ); }
}
/* info for indent_string()
func: indent a string including multiple lines
args: string str, int wrap-width, int firstline-indent, int extralines-indent
rtrn: string str
*/
string indent_string( string s_str,  int i_wrap, int i_indent1, int i_indent2 ){
string s_newstr, s_spaces;
if( i_wrap < 0 ){ i_wrap = 0; }
if( i_wrap > this_player()->GetScreen()[0] - 1 ){ i_wrap =  this_player()->GetScreen()[0] - 1; }
if( i_indent1 < 0 ){ i_indent1 = 0; }
if( i_indent2 < 0 ){ i_indent2 = 0; }
s_spaces =  "                                                                                ";
s_spaces += "                                                                                ";
s_spaces += "                                                                                ";
s_str = i_indent1 > 0 ? s_spaces[0..i_indent1] + s_str : s_str;
s_newstr = replace_string(
wrap( s_str, i_wrap ),
"\n",
"\n" + ( i_indent2 > 0 ? s_spaces[0..i_indent2] : "" )
)[0..<i_indent2];
return s_newstr[0..<4];
}
/* EOF */
« Last Edit: July 23, 2011, 01:06:51 AM by z993126 »

Offline z993126

  • BFF
  • ***
  • Posts: 128
    • View Profile
Re: Sky map and astronomical functions
« Reply #1 on: July 23, 2011, 01:10:53 AM »
Yet more files.
Code: [Select]
/* ** HEADER for ../funs/astro/chronfuns.c **
2005-Apr-19 - T. Cook
note: various functions from the following sources:
xylem.f2s.com/kepler
usno.gov
www.fourmilab.ch/documents/calendar
!! bounds for unix time 1901-Dec-13 to 2038-Jan-19 03:14:07 !!
-----
"type" (where relevant) argument values
0 = Gregorian "DDD, yyyy-MMM-dd HH:mm:ss UTC"   1 = Hebrew (T-="DDD, ")
2 = Islamic                                     3 = Persian
4 = Mayan (T="b.k.t.u.i, H, T")                 5 = Baha/'i (T+=" (k, v)"
6 = Indian Civil                                7 = Fr. Repub. (T="a-m-d-j")
8 = ISO8601a (T="yyyy-ww-d")                    9 = ISO8601b (T="yyyy-ddd")
10 = UNIX time()
-----
2011-Jul-17 - T. Cook began altering for DeadSouls
2011-Jul-13 - T. Cook fixed UTC2Local() to accomodate extended country code "%2s-%2s"
2010-Jul-05 - T. Cook changed UTC forms to "DDD, yyyy-MM-dd HH:mm:ss"
                     changed double-hypens to single
2010-Jun-11 - T. Cook corrected JD2UTC()
2010-Jun-10 - T. Cook added correction for unit rollovers in JD2UTC()
                     added monthDays()
                     added timezone corrections to DST()
                     modified DayOfYear() to use monthDays()
2010-Jun-02 - T. Cook added exception for three-letter Gregorian months in monthNumber()
2010-May-31 - T. Cook reformatted for notepad++
2010-Jan-24 - T. Cook fixed UTC2JD()
2010-Jan-12 - T. Cook tried to fix JD math
2009-Dec-30 - T. Cook removed unneeded JD2Gregorian functions
2009-Dec-20 - T. Cook started completely reworking to be general calendar converter
2009-Apr-15 - T. Cook began listing function dependencies & simplifying
2009-Mar-15 - T. Cook modified some
2009-Feb-11 - T. Cook begain modifying for tublib
2006-Apr-24 - T. Cook bug hunting
2005-Oct-30 - T. Cook modified DST() per 2007 changes
2005-Sep-28 - T. Cook tidied up
to do:
verify sidereal time functions
verify UTC to local Gregorian converter
add conversions for other calendars to/from JD
add day-of-year counting for other calendars
add formatting for crier data/time for other calendars
verify conversions at notable transitions (day/year/daylight-saving time)
add data for seasons for different latitudes/languages
add other countries' DST data
arrange functions by type (sub-day vs. regional calendar vs. astronomical/universal)
modify UTC2JD to allow pre-Gregorian dates
add data for locale Gregorian calendar change
add ability to use different-length days/subdivisions; make TimeOfDay() refer to sunpos()?
*/
#pragma strong_types

//------------------------------------------------------------------//
#include "/realms/ardneh/area/customdefs.h"
#include "/realms/ardneh/inc/astro.h"
#include "../../inc/math.h"

inherit MY_FUNS "/extras.c";
inherit MY_FUNS "/math.c";
inherit MY_FUNS "/lang/longnumbers.c";
inherit MY_FUNS "/geo/timezone.c";
//------------------------------------------------------------------//
string systime( mixed m_timezone, string s_geocode );
string Utime();
string UTC2Local( string s_UTC, mixed m_timezone, string s_geocode );
string JD2UTC( int i_JD, float f_dayfrac );
string JD2Local( int i_JD, float f_dayfrac, mixed m_timezone, string s_geocode );
mixed *Utime2JD();
string ampm( int i_hour );
float S( float f_T );
mixed *calc_seasons( int i_year );
mixed *season( int i_JD, float f_dayfrac, float f_latitude );
string TimeOfDay( int i_hour );
string MM2month( int i_month, int i_type );
int dayNumber( string s_day, int i_type );
int monthNumber( string s_month, int i_type );
int monthDays( int i_year, int i_month, int i_type );
int DayOfYear( string s_time, int i_type );
float mm2DayFrac( float f_minutes );
float DayFrac2mm( float f_dayfrac );
float DST( string s_time, float f_timezone, string s_geocode );
mixed *UTC2JD( string s_time );
string HH2HHmmss( float f_hours );
float mm2HH( float f_minutes );
float calcTimeJulianCent( float f_JulianDay );
float calcJDFromJulianCent( float f_t );
float JD2UTCmm( int i_JD, float f_dayfrac );
string mm2HHmmss( float f_minutes );
float HHmmss2mm( string s_time );
string HHmmss2hhmmss( string s_time );
string JD2GMST( int i_JD, float f_dayfrac );
string JD2LST( int i_JD, float f_dayfrac, float f_longitude );
string crierdate( string s_time, int i_type );
string criertime( string s_time );
/********************************************************************/
string *days = ({
({ // Julian & Gregorian day names
"Sunday", "Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday"
}),
({ // Hebrew day names
"Yom Rishon", "Yom Sheni", "Yom Shlishi", "Yom Revi'i", "Yom Chahmishi", "Yom Shishi", "Yom Shabbat"
}),
({ // Islamic day names
"al-'ahad", "al-'ithnayn", "ath-thalatha'", "al-'arb`a'", "al-khamis", "al-jum`a", "as-sabt"
}),
({ // Persian day names
"Yekshanbeh", "Doshanbeh", "Seshhanbeh", "Chaharshanbeh", "Panjshanbeh", "Jomeh", "Shanbeh"
}),
({ // Mayan day names
}),
({ // Baha/'i day names
"Jama//l", "Kama//l", "Fida//l", "Ida//l", "Istijla//l", "Istiqla//l", "Jala//l"
}),
({ // Indian Civil day names
"ravivara", "somavara", "mangalavara", "budhavara", "brahaspativara", "sukravara", "sanivara"
}),
({ // French Republic day names
}),
({ // ISO8601a day names
}),
({ // ISO8601b day names
}),
({ // UNIX time day names
})
});
string *months = ({
({ // Julian & Gregorian month names
"January", "February", "March", "April", "May", "June",
"July", "August", "September", "October", "November", "December"
}),
({ // Hebrew month names
"Nisan", "Iyyar", "Sivan", "Tammuz", "Av", "Elul",
"Tishri", "Heshvan", "Kislev", "Teveth", "Shevat", "Adar",
"Veadar"
}),
({ // Islamic month names
"Muharram", "Safar", "Rabi`al-Awwal", "Rabi`ath-Thani", "Jumada l-Ula", "Jumada t-Tania",
"Rajab", "Sha`ban", "Ramadan", "Shawwal", "Dhu l-Qa`da", "Dhu l-Hijja"
}),
({ // Persian month names
"Farvardin", "Ordibehesht", "Khordad", "Tir", "Mordad", "Shahrivar",
"Mehr", "Aban", "Azar", "Dey", "Bahman", "Esfand"
}),
({ // Mayan month names
}),
({ // Baha/'i month names
"Baha//", "Jala//l", "Jama//l", "`Azamat", "Nu//r", "Rahmat",
"Kalima//t", "Kala/l", "Asma//", "`Izzat", "Mashi/yyat", "`Ilm",
"Qudrat", "Qawl", "Masa//il", "Sharaf", "Sulta//n", "Mulk",
"Ayya//m-i-Ha//", "`Ala//'"
}),
({ // Indian Civil month names
"Caitra", "Vaisakha", "Jyaistha", "Asadha", "Sravana", "Bhadra",
"Asvina", "Kartika", "Agrahayana", "Pausa", "Magha", "Phalguna"
}),
({ // French Republic month names
"Vende/miaire", "Brumaire", "Frimaire", "Nivo^se", "Pluvio^se", "Vento^se",
"Germinal", "Flore/al", "Prairial", "Messidor", "Thermidor", "Fructidor",
"(Sans-culottides)"
}),
({ // ISO8601a month names
}),
({ // ISO8601b month names
}),
({ // UNIX time month names
})
});
string *years = ({
({ // Julian & Gregorian year names
}),
({ // Hebrew year names
}),
({ // Islamic year names
}),
({ // Persian year names
}),
({ // Mayan year names
}),
({ // Baha/'i year names
"Alif", "Ba//", "Ab", "Da//l", "Ba//b", "Va//v",
"Abad", "Ja//d", "Baha//", "Hubb", "Bahha//j", "Java//b",
"Ahad", "Vahha//b", "Vida//d", "Badi//", "Bahi", "Abha//",
"Vahi//d"
}),
({ // Indian Civil year names
}),
({ // French Republic year names
}),
({ // ISO8601a year names
}),
({ // ISO8601b year names
}),
({ // UNIX time year names
})
});
/* info for systime()
func: reformat the current system time value (Gregorian)
uses: ctime(), DST(), timezone()
args: float timezone offset, string country code
rtrn: string "DDD, yyyy-MM-dd HH:mm:ss TZ"
*/
string systime( mixed m_timezone, string s_geocode ){
string s_T, s_D, s_M, s_t;
float f_d;
s_T = ctime();
s_D = s_T[0..2];
s_M = s_T[4..6];
s_t = sprintf( "%s, %s-%s-%s %s", s_D, s_T[20..23], s_M, s_T[8..9], s_T[11..18] );
f_d = DST( s_t + " UTC", to_float( m_timezone ), s_geocode );
return sprintf( "%s %s", s_t, timezone( m_timezone, s_geocode, f_d )[0] );
}
/* info for Utime()
func: get the UTC time value (Gregorian)
uses: gmtime(), MM2month()
args: none
rtrn: string "DDD, yyyy-MMM-dd HH:mm:ss UTC"
*/
string Utime(){
mixed *ma_JD = Utime2JD();
return JD2UTC( ma_JD[0], ma_JD[1] );
}
/* info for UTC2Local()
func: convert UTC time value to local time Gregorian
('cheats' and converts to JD, offsets, converts to 'UTC')
args: string "DDD, yyyy-MMM-dd HH:mm:ss UTC
mixed zone offset, string country code
rtrn: string "DDD, yyyy-MMM-dd HH:mm:ss TZ"
*/
string UTC2Local( string s_UTC, mixed m_timezone, string s_geocode ){
float f_d = DST( s_UTC, to_float( m_timezone ), s_geocode[0..1] ), f_d2 = ( f_d + m_timezone ) / 24.0;
mixed *m_t = UTC2JD( s_UTC );
m_t[1] += f_d2;
if( m_t[1] >= 1 ){ ++m_t[0]; --m_t[1];}
if( m_t[1] <  0 ){ --m_t[0]; ++m_t[1];}
return JD2UTC( m_t[0], m_t[1] )[0..<4] + timezone( m_timezone, s_geocode[0..1], f_d )[0];
}
/* info for *Utime2JD()
func: Julian day plus fraction from UNIX time (cheaper)
uses: time()
args: none
rtrn: mixed array ({Julian day, day fraction})
*/
mixed *Utime2JD(){
int i_t = time();
float f_b = 0.5 + i_t % 86400 / 86400.0;
int i_a = to_int( i_t / 86400.0 ) + 2440586 + ( f_b < 0.5 ? 0 : 1 );
if( f_b >= 1 ){ --f_b; ++i_a; }
return ({ i_a, f_b });
}
/* info for JD2UTC()
func: convert Julian day plus fraction to UTC
uses: iwrap()
args: int Julian day
float day fraction
rtrn: string "DDD, yyyy-MMM-dd HH:mm:ss UTC"
*/
string JD2UTC( int i_JD, float f_dayfrac ){
float f_J = i_JD + f_dayfrac + 0.5;
int
i_j = to_int( f_J ) + 32044,
i_g = i_j / 146097,
i_dg = i_j % 146097,
i_c = ( i_dg / 36524 + 1 ) * 3 / 4,
i_dc = i_dg - i_c * 36524,
i_b = i_dc / 1461,
i_db = i_dc % 1461,
i_a = ( i_db / 365 + 1 ) * 3 / 4,
i_da = i_db - i_a * 365,
i_y = i_g * 400 + i_c * 100 + i_b * 4 + i_a,
i_m = ( i_da * 5 + 308 ) / 153 - 2,
i_d = i_da - ( i_m + 4 ) * 153 / 5 + 122,
i_Y = i_y - 4800 + ( i_m + 2 ) / 12,
i_M = ( ( i_m + 2 ) % 12 ) + 1,
i_D = i_d + 1;
float f_df = fwrap( ( f_dayfrac - 0.5 ) * 24.0, 0.0, 24.0 );
int
i_hr = iwrap( to_int( f_df ), 0, 24 ),
i_min = iwrap( to_int( ( f_df - i_hr ) * 60.0 ), 0, 60 ),
i_sec = iwrap( to_int( ( ( f_df - i_hr ) * 60.0 - i_min ) * 60.01 ), 0, 60 );
if( i_sec > 59 ){ i_sec = 0; ++i_min; }
if( i_min > 59 ){ i_min = 0; ++i_hr; }
if( i_hr > 23 ){ i_hr = 0; ++i_D; }
if( i_D > monthDays( i_Y, i_M, 0 ) ){ i_D = 1; ++i_M; }
if( i_M > 12 ){ i_M = 1; ++i_Y; }
return sprintf( "%s, %d-%s-%02d %02d:%02d:%02d UTC",
to_string( days[0][to_int( i_JD + f_dayfrac + 1.5 ) % 7] )[0..2],
i_Y,
to_string( months[0][i_M - 1])[0..2],
i_D,
i_hr,
i_min,
i_sec
);
}
/* info for JD2Local()
func: convert Julian Date value to local time Gregorian
args: int Julian Date, float day fraction, mixed zone offset, string country code
rtrn: string "DDD, yyyy-MMM-dd HH:mm:ss TZ"
*/
string JD2Local( int i_JD, float f_dayfrac, mixed m_timezone, string s_geocode ){
string s_UTC = JD2UTC( i_JD, f_dayfrac );
float f_d = DST( s_UTC, to_float( m_timezone ), s_geocode[0..1] ), f_d2 = ( f_d + m_timezone ) / 24.0;
mixed *m_t = ({ i_JD, f_dayfrac });//UTC2JD( s_UTC );
m_t[1] += f_d2;
if( m_t[1] >= 1 ){ ++m_t[0]; --m_t[1];}
if( m_t[1] <  0 ){ --m_t[0]; ++m_t[1];}
return JD2UTC( m_t[0], m_t[1] )[0..<4] + timezone( m_timezone, s_geocode[0..1], f_d )[0];
}
/* info for ampm()
func: get AM or PM status
uses: iwrap()
args: int hour
rtrn: string "AM", "PM"
*/
string ampm( int i_hour ){
return iwrap( i_hour, 0, 23 ) < 12 ? "AM" : "PM";
}
/* info for S()
func: maths for calc_seasons() function
uses: deg2rad()
args: float T (Julian century...?)
rtrn: some float
*/
float S( float f_T ){
return
485.0 * cos( deg2rad( 324.96 +   1934.136 * f_T ) ) +
203.0 * cos( deg2rad( 337.23 +  32964.467 * f_T ) ) +
199.0 * cos( deg2rad( 342.08 +     20.186 * f_T ) ) +
182.0 * cos( deg2rad(  27.85 + 445267.112 * f_T ) ) +
156.0 * cos( deg2rad(  73.14 +  45036.886 * f_T ) ) +
136.0 * cos( deg2rad( 171.52 +  22518.443 * f_T ) ) +
77.0 * cos( deg2rad( 222.54 +  65928.934 * f_T ) ) +
74.0 * cos( deg2rad( 296.72 +   3034.906 * f_T ) ) +
70.0 * cos( deg2rad( 243.58 +   9037.513 * f_T ) ) +
58.0 * cos( deg2rad( 119.81 +  33718.147 * f_T ) ) +
52.0 * cos( deg2rad( 297.17 +    150.678 * f_T ) ) +
50.0 * cos( deg2rad(  21.02 +   2281.226 * f_T ) ) +
45.0 * cos( deg2rad( 247.54 +  29929.562 * f_T ) ) +
44.0 * cos( deg2rad( 325.15 +  31555.956 * f_T ) ) +
29.0 * cos( deg2rad(  60.93 +   4443.417 * f_T ) ) +
18.0 * cos( deg2rad( 155.12 +  67555.328 * f_T ) ) +
17.0 * cos( deg2rad( 288.79 +   4562.452 * f_T ) ) +
16.0 * cos( deg2rad( 198.04 +  62894.029 * f_T ) ) +
14.0 * cos( deg2rad( 199.76 +  31436.921 * f_T ) ) +
12.0 * cos( deg2rad(  95.39 +  14577.848 * f_T ) ) +
12.0 * cos( deg2rad( 287.11 +  31931.756 * f_T ) ) +
12.0 * cos( deg2rad( 320.81 +  34777.259 * f_T ) ) +
  9.0 * cos( deg2rad( 227.73 +   1222.114 * f_T ) ) +
  8.0 * cos( deg2rad(  15.45 +  16859.074 * f_T ) );
}
/* info for *calc_seasons()
func: calculate season (equinox/solstice/midseason) Julian dates
uses: deg2rad(), S()
args: int year
rtrn: mixed array({
last year's equinox,
midseason,
last year's southern solstice
midseason,
equinox,
midseason,
northern solstice,
midseason,
equinox,
midseason,
southern solstice,
midseason
next year's equinox})
({date,day fraction})
note: converted to LPC from JS functions by Matt Oltersdorf
*/
mixed *calc_seasons( int i_year ){
int i_count;
float f_a, f_x, f_y, f_z, f_JDE0, *fa_s1, *fa_s2, f_T, f_W, f_dL;
mixed *ma_sT;
ma_sT = allocate( 13 );
for( i_count = 0; i_count < 13; ++i_count ){
ma_sT[i_count] = allocate( 2 );
}
f_a = to_float( i_year );
//write(sprintf("f_a=%O",f_a));
f_x = ( f_a - 1999.0 ) / 1000.0;
f_y = ( f_a - 2000.0 ) / 1000.0;
f_z = ( f_a - 2001.0 ) / 1000.0;
fa_s1 = allocate( 13 );
fa_s2 = allocate( 13 );
fa_s1[0]  = 2451810.21715 + ( 365242.01767 + ( 0.11575 - ( 0.00337 - 0.00078 * f_z ) * f_z ) * f_z ) * f_z;
fa_s1[2]  = 2451900.05952 + ( 365242.74049 + ( 0.06223 - ( 0.00823 - 0.00032 * f_z ) * f_z ) * f_z ) * f_z;
fa_s1[4]  = 2451623.80984 + ( 365242.37404 + ( 0.05169 - ( 0.00411 - 0.00057 * f_y ) * f_y ) * f_y ) * f_y;
fa_s1[6]  = 2451716.56767 + ( 365241.62603 + ( 0.00325 - ( 0.00888 - 0.00030 * f_y ) * f_y ) * f_y ) * f_y;
fa_s1[8]  = 2451810.21715 + ( 365242.01767 + ( 0.11575 - ( 0.00337 - 0.00078 * f_y ) * f_y ) * f_y ) * f_y;
    fa_s1[10] = 2451900.05952 + ( 365242.74049 + ( 0.06223 - ( 0.00823 - 0.00032 * f_y ) * f_y ) * f_y ) * f_y;
fa_s1[12] = 2451623.80984 + ( 365242.37404 + ( 0.05169 - ( 0.00411 - 0.00057 * f_x ) * f_x ) * f_x ) * f_x;
for( i_count = 0; i_count < 13; i_count += 2 ){ // generate midseason times
//write(sprintf("fa_s1[%d]=%O",i_count,fa_s1[i_count]));
f_T = ( fa_s1[i_count] - 2451545.0 ) / 36525.0;
//write(sprintf("f_T[%d]=%O",i_count,f_T));
f_W = deg2rad( 35999.373 * f_T - 2.47 );
//write(sprintf("f_W[%d]=%O",i_count,f_W));
f_dL = 1.0 + 0.0334 * cos( f_W ) + 0.0007 * cos( f_W );
//write(sprintf("f_dL[%d]=%O",i_count,f_dL));
fa_s2[i_count] = fa_s1[i_count] + ( 0.00001 * S( f_T ) ) / f_dL - ( 66.0 + ( i_year - 2000 ) * 1.0 ) / 86400.0;
//write(sprintf("fa_s2[%d]=%O",i_count,fa_s2[i_count]));
ma_sT[i_count][0] = to_int( fa_s2[i_count] ); // Julian date of event
ma_sT[i_count][1] = fa_s2[i_count] - ma_sT[i_count][0]; // day fraction of event
//write("-----");
}
for( i_count = 1; i_count < 13; i_count += 2 ){ // generate season-change times
ma_sT[i_count][0] = to_int( ( fa_s2[i_count - 1] + fa_s2[i_count + 1] ) / 2.0 );
ma_sT[i_count][1] = ( ( fa_s2[i_count - 1] + fa_s2[i_count + 1] ) / 2.0 ) - ma_sT[i_count][0];
}
return ma_sT;
}
/* info for *season()
func: get the season name at UTC for given Julian date (English)
args: int Julian date
float day fraction
float latitude
rtrn: mixed array ({ string season, float began, float ends })
*/
mixed *season( int i_JD, float f_dayfrac, float f_latitude ){
float f_d = i_JD + f_dayfrac;
string s_b, s_c, s_name;
int i_count;
int i_y;
mixed *ma_a, *ma_s;
float
f_TROPICN = 23.4381,
f_TROPICS = -f_TROPICN,
f_POLARN = 90.0 - f_TROPICN,
f_POLARS = -f_POLARN;
sscanf( JD2UTC( i_JD, f_dayfrac ), "%s, %d-%s UTC", s_b, i_y, s_c );
ma_a = calc_seasons( i_y );
ma_s = allocate( 13 );
for( i_count = 0; i_count < 13; ++i_count ){ ma_s[i_count] = ma_a[i_count][0] + ma_a[i_count][1]; }
if( f_d >= ma_s[3] && f_d < ma_s[5] ){
if( f_latitude >= f_POLARN ){                                s_name = "twilight"; }
else if( f_latitude >= f_TROPICN && f_latitude < f_POLARN ){ s_name = "spring"; }
else if( f_latitude >= 0 && f_latitude < f_TROPICN ){        s_name = "wet"; }
else if( f_latitude < 0 && f_latitude >= f_TROPICS ){        s_name = "dry"; }
else if( f_latitude < f_TROPICS && f_latitude >= f_POLARS ){ s_name = "autumn"; }
else{                                                        s_name = "twilight"; }
return ({ s_name, ma_s[3], ma_s[5] });
}else if( f_d >= ma_s[5] && f_d < ma_s[7] ){
if( f_latitude >= f_POLARN ){                                s_name = "midnight sun"; }
else if( f_latitude >= f_TROPICN && f_latitude < f_POLARN ){ s_name = "summer"; }
else if( f_latitude >= 0 && f_latitude < f_TROPICN ){        s_name = "wet";}
else if( f_latitude < 0 && f_latitude >= f_TROPICS ){        s_name = "dry"; }
else if( f_latitude < f_TROPICS && f_latitude >= f_POLARS ){ s_name = "winter"; }
else{                                                        s_name = "polar night"; }
return ({ s_name, ma_s[5], ma_s[7] });
}else if( f_d >= ma_s[7] && f_d < ma_s[9] ){
if( f_latitude >= f_POLARN ){                                s_name = "twilight"; }
else if( f_latitude >= f_TROPICN && f_latitude < f_POLARN ){ s_name = "autumn"; }
else if( f_latitude >= 0 && f_latitude < f_TROPICN ){        s_name = "dry"; }
else if( f_latitude < 0 && f_latitude >= f_TROPICS ){        s_name = "wet"; }
else if( f_latitude < f_TROPICS && f_latitude >= f_POLARS ){ s_name = "spring"; }
else{                                                        s_name = "twilight"; }
return ({ s_name, ma_s[7], ma_s[9] });
}else{
if( f_latitude >= f_POLARN ){                                s_name = "polar night"; }
else if( f_latitude >= f_TROPICN && f_latitude < f_POLARN ){ s_name = "winter"; }
else if( f_latitude >= 0 && f_latitude < f_TROPICN ){        s_name = "dry"; }
else if( f_latitude < 0 && f_latitude >= f_TROPICS ){        s_name = "wet"; }
else if( f_latitude < f_TROPICS && f_latitude >= f_POLARS ){ s_name = "summer"; }
else{                                                        s_name = "midnight sun"; }
if( f_d >= ma_s[9] ){ return ({ s_name, ma_s[9], ma_s[11] }); }
else{ return ({ s_name, ma_s[1], ma_s[3] }); }
}
}
/* info for TimeOfDay()
func: get the general time of 24-hour day (English)
args: int hour
rtrn: string "early morning", "late morning", "afternoon", "evening", "tempus incognito"
*/
string TimeOfDay( int i_hour ){
switch( iwrap( i_hour, 0, 23 )){
case 00..05: return "early morning";
case 06..11: return "late morning";
case 12..17: return "afternoon";
case 18..23: return "evening";
default:     return "tempus incognito";
}
}
/* info for MM2month()
func: get the short month name from the month number
args: int month, int type (see chart @ top)
rtrn: string shortform month
*/
string MM2month( int i_month, int i_type ){
return to_string( months[i_type][i_month - 1] )[0..2];
}
/* info for monthNumber()
func: get the number of the month
args: string month, int type (see table @ top)
rtrn: int month number
*/
int monthNumber( string s_month, int i_type ){
int i_count, i_s;
for(i_count = 0, i_s = sizeof( months[i_type] ); i_count < i_s; ++i_count ){
if( i_type == 0 ){
if( lower_case( s_month[0..2] ) == lower_case( to_string( months[i_type][i_count] ) )[0..2] ){
return ++i_count;
}
}else if( lower_case( s_month ) == lower_case( to_string( months[i_type][i_count] ) ) ){
return ++i_count;
}
}
return 1;
}
/* info for dayNumber()
func: get the number of the weekday (Sunday=1)
args: string day, int type (see table @ top)
rtrn: int weekday number
*/
int dayNumber( string s_day, int i_type ){
int i_count, i_s;
for( i_count = 0, i_s = sizeof( days[i_type] ); i_count < i_s; ++i_count ){
if( i_type == 0 ){
if( lower_case( s_day[0..2] ) == lower_case( to_string( days[i_type][i_count] ) )[0..2] ){
return ++i_count;
}
}else if( lower_case( s_day ) == lower_case( to_string( days[i_type][i_count] ) ) ){
return ++i_count;
}
}
return 1;
}
/* info for monthDays()
func: return number of days in given month (only Gregorian for now)
args: int year, int month number, type (see table @ top)
rtrn: int number of days in month [for Gregorian]
*/
int monthDays( int i_year, int i_month, int i_type ){
switch( i_type ){
case 0:
switch( i_month ){
case 9: case 4: case 6: case 11: return 30;
case 2:
if( ( i_year % 4 == 0 && i_year % 100 !=0 ) || i_year % 400 == 0 ){ return 29; }else{ return 28; }
break;
case 1: case 3: case 5: case 7: case 8: case 10: case 12: return 31;
default: return 0;
}
default: return 0;
}
}
/* info for DayOfYear()
func: calculate the day of year (only Gregorian for now)
args: string "yyyy-MMM-dd", type (see table @ top)
rtrn: int day of year 1 .. 365 (or 366) [for Gregorian]
*/
int DayOfYear( string s_time, int i_type ){
int i_y, i_d, i_count, i_m, i_D;
string s_M;
switch( i_type ){
case 0:
if( sscanf( s_time, "%d-%s-%([0-9][0-9])", i_y, s_M, i_d ) != 3 ){
write( "Error in input to DayOfYear() - '" + s_time + "'.\n" );
return 0;
}else{
i_d = to_int( i_d );
i_m = monthNumber( s_M, i_type ), i_D = i_d;
for( i_count = 1; i_count < 12; ++i_count ){
if( i_m > i_count ){ i_D += monthDays( i_y, i_count, 0 ); }
}
return i_D;
}
default: return 0;
}
}
/* info for mm2DayFrac()
func: convert time in minutes to fraction of day
args: float time in minutes
rtrn: float fraction of day
*/
float mm2DayFrac( float f_minutes ){
return to_float( f_minutes ) / 1440.0;
}
/* info for DayFrac2mm()
func: convert fraction of day to time in minutes
args: float fraction of day
rtrn: float time in minutes
*/
float DayFrac2mm( float f_dayfrac ){
return f_dayfrac * 1440.0;
}
/* info for DST()
func: calculate Daylight Saving Time status from Gregorian date UTC
args: string "DDD, yyyy-MMM-dd HH:mm:ss UTC"
float time zone
string country code
rtrn: float DST offset in hours
*/
float DST( string s_time, float f_timezone, string s_geocode ){
int i_y, i_d, i_H, i_m, i_s, i_DSTStart, i_DSTEnd, i_a1, i_month;
string s_D, s_M;
if( sscanf( s_time, "%s, %d-%s-%([0-9][0-9]) %([0-9][0-9]):%([0-9][0-9]):%([0-9][0-9]) UTC", s_D, i_y, s_M, i_d, i_H, i_m, i_s ) != 7 ){
write( "Error in input to DST() - '" + s_time + "'.\n" );
return 0.0;
}else{
i_d = to_int(i_d);
i_H = to_int(i_H);
i_m = to_int(i_m);
i_s = to_int(i_s);
i_month = monthNumber( s_M, 0 );
// begin fiddly bit to correct for time zone difference
if( f_timezone - to_int( f_timezone ) != 0 ){
i_m += to_int( 60 * ( f_timezone - to_int( f_timezone ) ) );
if( i_m < 0 ){ --i_H; }
if( i_m > 59 ){ ++i_H; }
i_m = iwrap( i_m, 0, 59 );
}
i_H += to_int( f_timezone );
if( i_H < 0 || i_H > 23 ){
if( i_H < 0 ){ --i_d; }
if( i_H > 23 ){ ++i_d; }
i_H = iwrap( i_H, 0, 23 );
}
if( i_d < 1 || i_d > monthDays( i_y, i_month, 0 ) ){
if( i_d < 1 ){
--i_month;
}
if( i_d > monthDays( i_y, i_month, 0 ) ){
i_d -= monthDays( i_y, i_month, 0 ); ++i_month;
}
i_month = iwrap( i_month, 1, 12 );
i_d = iwrap( i_d, 1, monthDays( i_y, i_month, 0 ) );
}
if( i_month < 1 || i_month > 12 ){
if( i_month < 1 ){ --i_y; }
if( i_month > 12){ ++i_y; }
i_month = iwrap( i_month, 1, 12 );
}
// end time zone correction
switch( lower_case( s_geocode[0..1] ) ){
case "us":
switch( i_y ){
case 1976..2006:
switch( i_month ){
case 4: // first Sunday of Apr.
i_a1 = ( 2 + 6 * i_y - ( i_y / 4 ) ) % 7 + 1;
if( i_d == i_a1 ){ if( i_H >= 2 ){ return 1.0; }else{ return 0.0; } }
else if( i_d > i_a1 ){ return 1.0; }
else{ return 0.0; }
case 5..9: return 1.0;
case 10: // last Sunday of Oct.
i_a1 = ( 31 - ( ( i_y * 5 / 4 ) + 1 ) % 7 );
if( i_d == i_a1 ){ if( i_H >= 1 ){ return 0.0; }else{ return 1.0; } }
else if( i_d > i_a1 ){ return 0.0; }
else{ return 1.0; }
default: return 0.0;
}
case 2007..2099:
switch( i_month ){
case 3: // second Sunday of Mar.
i_a1 = ( ( 5 + 6 * i_y - ( i_y / 4 ) ) % 7 + 8 );
if( i_d == i_a1 ){ if( i_H >= 2 ){ return 1.0; }else{ return 0.0; } }
else if( i_d > i_a1 ){ return 1.0; }
else{ return 0.0; }
case 4..10: return 1.0;
case 11: // first Sunday of Nov.
i_a1 = ( 7 - ( ( i_y * 5 / 4 ) + 1 ) % 7 );
if( i_d == i_a1 ){ if( i_H >= 1 ){ return 0.0; }else{ return 1.0; } }
else if( i_d > i_a1 ){ return 0.0; }
else{ return 1.0; }
default: return 0.0;
}
default: return 0.0;
}
default: return 0.0;
}
}
}
/* info for *UTC2JD()
func: convert Gregorian date UTC to Julian day
args: string "DDD, yyyy-MMM-dd HH:mm:ss UTC"
rtrn: mixed array ({int Julian day, float day fraction})
note: valid for dates after 1752 Sep 01
*/
mixed *UTC2JD( string s_time ){
int i_y, i_d, i_H, i_mm, i_s, i_mon, i_a, i_yr, i_m, i_JDN;
float df;
string s_D, s_M;
mixed *ma_JDArray;
if( sscanf( s_time, "%s, %d-%s-%([0-9][0-9]) %([0-9][0-9]):%([0-9][0-9]):%([0-9][0-9]) UTC", s_D, i_y, s_M, i_d, i_H, i_mm, i_s ) != 7 ){
write( "Error in input to UTC2JD() - '" +  s_time + "'.\n" );
return ({ 0, 0.0 });
}else{
/* alt method:
int mmm=monthNumber(M,0);
int b=(1461*(y+4800+(mmm-14)/12))/4+
(367*(mmm-2-12*((mmm-14)/12)))/12-
(3*((y+4900+(mmm-14)/12)/100))/4+d-32076;
*/
ma_JDArray = ({ 0, 0.0 });
i_d = to_int( i_d );
i_H = to_int( i_H );
i_mm = to_int( i_mm );
i_s = to_int( i_s );
i_mon = monthNumber( s_M, 0 );
i_a = ( 14 - i_mon ) / 12;
i_yr = i_y + 4800 - i_a;
i_m = i_mon + ( 12 * i_a ) - 3;
i_JDN = i_d + ( ( ( 153 * i_m ) + 2 ) / 5 ) + ( 365 * i_yr ) + ( i_yr / 4 ) - ( i_yr / 100 ) + ( i_yr / 400 ) - 32045;
ma_JDArray[0] = to_int( i_JDN );
df = HHmmss2mm( sprintf( "%02d:%02d:%02d", iwrap(( i_H - 12 ),0,23), i_mm, i_s ) ) / 1440.0;
ma_JDArray[1] = df;
if( ma_JDArray[1] <  0.0 ){ --ma_JDArray[0]; ++ma_JDArray[1]; }
if( ma_JDArray[1] >= 0.5 ){ --ma_JDArray[0]; }
return ma_JDArray;
}
}
/* info for HH2HHmmss()
func: convert 24-hour plus fractional hour to string "HH:mm:ss"
uses: to_int(), sprintf()
args: float hour 0..<24
rtrn: string "HH:mm:ss"
*/
string HH2HHmmss( float f_hours ){
int i_H, i_m, i_s;
i_H = to_int( f_hours );
i_m = to_int( ( f_hours - i_H ) * 60.00 );
i_s = to_int( ( ( ( f_hours - i_H ) * 60.00 ) - i_m ) * 60.00 );
return sprintf( "%02d:%02d:%02d", i_H, i_m, i_s );
}
/* info for mm2HH
func: convert the time in minutes to 24-hour plus fraction
uses: to_float()
args: float time in minutes
rtrn: float 0..<24
*/
float mm2HH( float f_minutes ){
return f_minutes / 60.0;
}
/* info for calcTimeJulianCent()
func: calculate the Julian century since J2000.0
args: float Julian day
rtrn: float Julian century T since J2000.0
*/
float calcTimeJulianCent( float f_JulianDay ){
return ( f_JulianDay - 2451545.0 ) / 36525.0;
}
/* info for calcJDFromJulianCent()
func: convert centures since J2000.0 to Julian day
args: float number of Julian centuries since J2000.0
rtrn: float Julian day
*/
float calcJDFromJulianCent( float f_t ){
return f_t * 36525.0 + 2451545.0;
}
/* info for JD2UTCmm()
func: convert Julian date to UTC minutes
uses: calcHourFromJD()
args: int Julian day
float day fraction
rtrn: float UTC time in minutes
*/
float JD2UTCmm( int i_JD, float f_dayfrac ){
int i_H, i_m, i_s;
string s_a;
sscanf( JD2UTC( i_JD, f_dayfrac ), "%s %([0-9][0-9]):%([0-9][0-9]):%([0-9][0-9]) UTC", s_a, i_H, i_m, i_s );
i_H = to_int( i_H );
i_m = to_int( i_m );
i_s = to_int( i_s );
return ( i_H * 60.0 + i_m + i_s / 60.0 ) / 1440.0;
}
/* info for mm2HHmmss()
func: convert time in minutes to twenty-four hour time-code
args: float time in minutes
rtrn: string time "HH:mm:ss"
*/
string mm2HHmmss( float f_minutes ){
int i_H, i_m, i_s;
i_H = to_int( f_minutes / 60.00 );
i_m = to_int( f_minutes ) % 60;
i_s = to_int( ( f_minutes - to_int( f_minutes ) ) * 60.00 );
return sprintf( "%02d:%02d:%02d", i_H, i_m, i_s );
}
/* info for HHmmss2mm()
func: convert time code to time fraction in minutes
args: string T "HH:mm:ss"
rtrn: float time in minutes
*/
float HHmmss2mm( string s_time ){
int i_H, i_m, i_s;
float f_mm;
if( sscanf( s_time, "%d:%([0-9][0-9]):%([0-9][0-9])", i_H, i_m, i_s ) != 3 ){ return 0.0; }else{
i_m = to_int( i_m );
i_s = to_int( i_s );
f_mm = i_H * 60.00 + i_m + i_s / 60.00;
f_mm = fwrap( f_mm, 0.0, 1440.0 );
return f_mm;
}
}
/* info for HHmmss2hhmmss()
func: convert twenty-four to twelve-hour time
args: string T "HH:mm:ss"
rtrn: string T "hh:mm:ss ?M"
*/
string HHmmss2hhmmss( string s_time ){
int i_h, i_m, i_s;
if( sscanf( s_time, "%([0-9][0-9]):%([0-9][0-9]):%([0-9][0-9])", i_h, i_m, i_s ) != 3 ){
return "tempus igcognito";
}else{
i_h = to_int( i_h );
i_m = to_int( i_m );
i_s = to_int( i_s );
return sprintf( "%02d:%02d:%02d %2s", iwrap( i_h, 1, 12 ), i_m, i_s, ampm( i_h ) );
}
}
/* info for JD2GMST()
func: Greenwich Mean Sidereal Time
args: int JD, float dayfrac
rtrn: string "HH:mm:ss GMST"
*/
string JD2GMST( int i_JD, float f_dayfrac ){
string s_a, s_UM;
int i_Uy, i_Ud, i_UH, i_Um, i_Us, i_UMM;
float f_dj, f_T, f_Udec, f_R1, f_R0, f_t0, f_gstdec;
sscanf( JD2UTC( i_JD, f_dayfrac ), "%s, %d-%s-%([0-9][0-9]) %([0-9][0-9]):%([0-9][0-9]):%([0-9][0-9]) UTC", s_a, i_Uy, s_UM, i_Ud, i_UH, i_Um, i_Us );
i_Ud = to_int( i_Ud );
i_UH = to_int( i_UH );
i_Um = to_int( i_Um );
i_Us = to_int( i_Us );
i_UMM = monthNumber( s_UM, 0 );
f_dj = ( i_JD + f_dayfrac ) - 2415020;
f_T = ( f_dj / 36525.0 ) - 1.0;
f_Udec = ( ( i_Us / 60.0 ) + i_Um ) / 60.0 + i_UH;
f_R1 = 6.697374558 + ( 2400.0 * ( f_T - ( ( i_Uy - 2000 ) / 100 ) ) );
f_R0 = ( 0.0513366 * f_T ) + ( 0.00002586222 * f_T * f_T ) - ( 0.000000001722 * f_T * f_T * f_T );
f_t0 = fwrap( f_R0 + f_R1, 0.0, 24.0 );
f_gstdec = fwrap( ( f_Udec * 1.002737935 ) + f_t0, 0.0, 24.0 );
i_UH = to_int( f_gstdec );
i_Um = to_int( ( f_gstdec - i_UH ) * 60.0 );
i_Us = to_int( ( ( f_gstdec - i_UH ) * 60.0 - i_Um ) * 60.0 );
return sprintf( "%02d:%02d:%02d GMST", i_UH, i_Um, i_Us );
}
/* info for JD2LST()
func: calculate the local sidereal time
args: int Julian day
float dayfrac
float longitude in degrees (west negative)
rtrn: string sidereal time "HH:mm:ss LST"
*/
string JD2LST( int i_JD, float f_dayfrac, float f_longitude ){
string s_a, s_UM;
int i_Uy, i_Ud, i_UH, i_Um, i_Us, i_UMM;
float f_dj, f_T, f_Udec, f_R1, f_R0, f_t0, f_gstdec;
sscanf( JD2UTC( i_JD, f_dayfrac ), "%s, %d-%s-%([0-9][0-9]) %([0-9][0-9]):%([0-9][0-9]):%([0-9][0-9]) UTC", s_a, i_Uy, s_UM, i_Ud, i_UH, i_Um, i_Us );
i_Ud = to_int( i_Ud );
i_UH = to_int( i_UH );
i_Um = to_int( i_Um );
i_Us = to_int( i_Us );
i_UMM = monthNumber( s_UM, 0 );
f_dj = ( i_JD + f_dayfrac ) - 2415020;
f_T = ( f_dj / 36525.0 ) - 1.0;
f_Udec = ( ( i_Us / 60.0 ) + i_Um ) / 60.0 + i_UH;
f_R1 = 6.697374558 + ( 2400.0 * ( f_T - ( ( i_Uy - 2000 ) / 100 ) ) );
f_R0 = ( 0.0513366 * f_T ) + ( 0.00002586222 * f_T * f_T ) - ( 0.000000001722 * f_T * f_T * f_T );
f_t0 = fwrap( f_R0 + f_R1, 0.0, 24.0 );
f_gstdec = fwrap( ( f_Udec * 1.002737935 ) + f_t0 + ( f_longitude / 15.0 ), 0.0, 23.99999 );
i_UH = to_int( f_gstdec );
i_Um = to_int( ( f_gstdec - i_UH ) * 60.0 );
i_Us = to_int( ( ( f_gstdec - i_UH) * 60.0 - i_Um ) * 60.0 );
return sprintf( "%02d:%02d:%02d LST", i_UH, i_Um, i_Us );
}
/* info for crierdate()
func: format the date into words (English)
args: string "DDD, yyyy-MMM-dd HH:mm:ss", int type (see list @ top) [semi-ignored for now]
rtrn: string [day], the [n][st/nd/th] of [month], [year]
*/
string crierdate( string s_time, int i_type ){
string s_D, s_M, s_x;
int i_y, i_d;
if( sscanf( s_time, "%s, %d-%s-%([0-9][0-9]) %s", s_D, i_y, s_M, i_d, s_x ) != 5 ){
return "an unknown date in time";
}else{
i_d = to_int( i_d );
return sprintf(
"%s, the %s of %s, in the %s year %s the Common Era",
days[i_type][dayNumber( s_D, i_type ) - 1],
spell_number( i_d, 1 ),
months[i_type][monthNumber( s_M, i_type ) - 1],
spell_number( abs( i_y ), 1 ),
i_y < 0 ? "before" : "of"
);
}
}
/* info for criertime()
func: format the time into words (English)
args: string "hh:mm"
rtrn: string [minute(s) [past/'til]] [hour] o'clock in the [time of day]
*/
string criertime( string s_time ){
int i_H, i_h, i_m;
string s_longtime;
if( sscanf( s_time, "%([0-9][0-9]):%([0-9][0-9])", i_H, i_m ) != 2 ){ return "eternity"; }else{
i_H = to_int( i_H );
i_m = to_int( i_m );
i_h = iwrap( i_H, 1, 12 ); i_m = iwrap( i_m, 0, 59 );
s_longtime = "";
if( i_m < 31 ){
switch( i_m ){
case 0: break;
case 1: s_longtime += "a minute"; break;
case 15: s_longtime += "a quarter"; break;
case 30: s_longtime += "half"; break;
default: s_longtime += spell_number( i_m, 0 ) + " minutes"; break;
}
if( i_m != 0 ){ s_longtime += " past "; }
}else{
i_h = iwrap( ++i_h, 1, 12 );
switch( i_m ){
case 45: s_longtime += "a quarter"; break;
case 59: s_longtime + "a minute"; break;
default: s_longtime += spell_number( 60 - i_m, 0 ) + " minutes"; break;
}
s_longtime += " 'til ";
}
s_longtime += spell_number( i_h, 0 ) + " o'clock in the " + TimeOfDay( i_H );
return s_longtime;
}
}
/* EOF */

Offline z993126

  • BFF
  • ***
  • Posts: 128
    • View Profile
Re: Sky map and astronomical functions
« Reply #2 on: July 23, 2011, 01:21:22 AM »
Code: [Select]
/* ** HEADER for ../funs/boxdrawing.c **
2005-Jul-31 - T. Cook
boxdrawing character handling
-----
2011-Jul-17 - T. Cook began altering for DeadSouls
2010-May-21 - T. Cook reformatted for notepad++
2009-Dec-30 - T. Cook modified for tublib, streamlined
2005-Sep-27 - T. Cook tidied up
*/
#pragma strong_types

string box(int form,int style);
/* info for box()
func: get boxdrawing character
args: int form (UDRL  1=single, 2=double, 0=none)
int style (1=OEM, 2=UTF-8, 0=ASCII [default])
rtrn: string boxdrawing character
*/
string box(int form,int style){
switch( style ){
case 1:
switch( form ){
case 0000: return "\032";
case 0011: return "\196";
case 0022: return "\205";
case 0101: return "\191";
case 0102: return "\184";
case 0110: return "\218";
case 0111: return "\194";
case 0120: return "\213";
case 0122: return "\209";
case 0201: return "\183";
case 0202: return "\187";
case 0210: return "\214";
case 0211: return "\210";
case 0220: return "\201";
case 0222: return "\203";
case 1001: return "\217";
case 1002: return "\190";
case 1010: return "\192";
case 1011: return "\193";
case 1020: return "\212";
case 1022: return "\207";
case 1100: return "\179";
case 1101: return "\180";
case 1102: return "\181";
case 1110: return "\195";
case 1111: return "\197";
case 1120: return "\198";
case 1122: return "\216";
case 2001: return "\189";
case 2002: return "\188";
case 2010: return "\211";
case 2011: return "\208";
case 2020: return "\200";
case 2022: return "\202";
case 2200: return "\186";
case 2201: return "\182";
case 2202: return "\185";
case 2210: return "\199";
case 2211: return "\215";
case 2220: return "\204";
case 2222: return "\206";
default:   return "\063";
}
break;
case 2:
switch( form ){
case 0000: return "\227\128\128";
case 0011: return "\226\148\128";
case 0022: return "\226\149\144";
case 0101: return "\226\148\144";
case 0102: return "\226\149\149";
case 0110: return "\226\148\140";
case 0111: return "\226\148\172";
case 0120: return "\226\149\146";
case 0122: return "\226\149\164";
case 0201: return "\226\149\150";
case 0202: return "\226\149\151";
case 0210: return "\226\149\147";
case 0211: return "\226\149\165";
case 0220: return "\226\149\148";
case 0222: return "\226\149\166";
case 1001: return "\226\148\152";
case 1002: return "\226\149\155";
case 1010: return "\226\148\148";
case 1011: return "\226\148\180";
case 1020: return "\226\149\152";
case 1022: return "\226\149\167";
case 1100: return "\226\148\130";
case 1101: return "\226\148\164";
case 1102: return "\226\149\161";
case 1110: return "\226\148\156";
case 1111: return "\226\148\188";
case 1120: return "\226\149\158";
case 1122: return "\226\149\170";
case 2001: return "\226\149\156";
case 2002: return "\226\149\157";
case 2010: return "\226\149\153";
case 2011: return "\226\149\168";
case 2020: return "\226\149\154";
case 2022: return "\226\149\169";
case 2200: return "\226\149\145";
case 2201: return "\226\149\162";
case 2202: return "\226\149\163";
case 2210: return "\226\149\159";
case 2211: return "\226\149\171";
case 2220: return "\226\149\160";
case 2222: return "\226\149\172";
default:   return "\239\188\159";
}
break;
default:
switch( form ){
case 0000: return " ";
case 0011: return "-";
case 0022: return "=";
case 0101: return "+";
case 0102: return "+";
case 0110: return "+";
case 0111: return "+";
case 0120: return "+";
case 0122: return "+";
case 0201: return "+";
case 0202: return "+";
case 0210: return "+";
case 0211: return "+";
case 0220: return "+";
case 0222: return "+";
case 1001: return "+";
case 1002: return "+";
case 1010: return "+";
case 1011: return "+";
case 1020: return "+";
case 1022: return "+";
case 1100: return "|";
case 1101: return "+";
case 1102: return "+";
case 1110: return "+";
case 1111: return "+";
case 1120: return "+";
case 1122: return "+";
case 2001: return "+";
case 2002: return "+";
case 2010: return "+";
case 2011: return "+";
case 2020: return "+";
case 2022: return "+";
case 2200: return "|";
case 2201: return "+";
case 2202: return "+";
case 2210: return "+";
case 2211: return "+";
case 2220: return "+";
case 2222: return "+";
default:   return "?";
}
break;
}
return "";
}
/* EOF */

Code: [Select]
/* ** HEADER for ../funs/geo/regional_divisions.c **
2010-Jun-01 - T. Cook
get names of political divisions from abbreviations
-----
2011-Jul-17 - T. Cook began altering for DeadSouls
to do:
finish adding data to region_name
*/
#pragma strong_types

/* info for countryname()
func: expand country names
args: string 2-letter ISO 3166 country code
rtrn: string country name
*/
string countryname(string cc){
switch(lower_case(cc)){
case "ad":return "Andorra";
case "ac":return "Ascension Island (reserved)";
case "ae":return "the United Arab Emirates";
case "af":return "Afghanistan";
case "ag":return "Antigua and Barbuda";
case "ai":return "Anguilla";
case "al":return "Albania";
case "am":return "Armenia";
case "an":return "the Netherlands Antilles";
case "ao":return "Angola";
case "aq":return "Antarctica";
case "ar":return "Argentina";
case "as":return "American Samoa";
case "at":return "Austria";
case "au":return "Australia";
case "aw":return "Aruba";
case "ax":return "A^land Islands";
case "az":return "Azerbaijan";

case "ba":return "Bosnia and Herzegowina";
case "bb":return "Barbados";
case "bd":return "Bangladesh";
case "be":return "Belgium";
case "bf":return "Burkina Faso";
case "bg":return "Bulgaria";
case "bh":return "Bahrain";
case "bi":return "Burundi";
case "bj":return "Benin";
case "bl":return "Saint Barthe/lemy";
case "bm":return "Bermuda";
case "bn":return "Brunei Darussalam";
case "bo":return "the Plurinational State of Bolivia";
case "br":return "Brazil";
case "bs":return "the Bahamas";
case "bt":return "Bhutan";
case "bu":return "Burma (obsolete)";
case "bv":return "Bouvet Island";
case "bw":return "Botswana";
case "by":return "Belarus";
case "bz":return "Belize";

case "ca":return "Canada";
case "cc":return "the Cocos (Keeling) Islands";
case "cd":return "the Democratic Republic of Congo (was Zaire)";
case "cf":return "the Central African Republic";
case "cg":return "the People's Republic of Congo";
case "ch":return "Switzerland";
case "ci":return "Co^te d'Ivoire";
case "ck":return "the Cook Islands";
case "cl":return "Chile";
case "cm":return "Cameroon";
case "cn":return "People's Republic of China";
case "co":return "Colombia";
case "cp":return "Clipperton Island (reserved)";
case "cr":return "Costa Rica";
case "cs":return "Czechoslovakia or Serbia and Montenegro (obsolete)";
case "cu":return "Cuba";
case "cv":return "Cape Verde";
case "cx":return "Christmas Island";
case "cy":return "Cyprus";
case "cz":return "the Czech Republic";

case "de":return "Germany";
case "dg":return "Diego Garcia (reserved)";
case "dj":return "Djibouti";
case "dk":return "Denmark";
case "dm":return "Dominica";
case "do":return "the Dominican Republic";
case "dz":return "Algeria";

case "ea":return "Melilla Ceuta (reserved)";
case "ec":return "Ecuador";
case "ee":return "Estonia";
case "eg":return "Egypt";
case "eh":return "Western Sahara";
case "er":return "Eritrea";
case "es":return "Spain";
case "et":return "Ethiopia";
case "eu":return "the European Union (reserved)";

case "fi":return "Finland";
case "fj":return "Fiji";
case "fk":return "the Falkland Islands (Malvinas)";
case "fm":return "the Federated States of Micronesia";
case "fo":return "the Faroe Islands";
case "fr":return "France";
case "fx":return "Metropolitan France (reserve)";

case "ga":return "Gabon";
case "gb":return "the United Kingdom";
case "gd":return "Grenada";
case "ge":return "Georgia";
case "gf":return "French Guiana";
case "gg":return "Guernsey";
case "gh":return "Ghana";
case "gi":return "Gibraltar";
case "gl":return "Greenland";
case "gm":return "the Republic of the Gambia";
case "gn":return "Guinea";
case "gp":return "Guadeloupe";
case "gq":return "Equatorial Guinea";
case "gr":return "Greece";
case "gs":return "South Georgia and the South Sandwich Islands";
case "gt":return "Guatemala";
case "gu":return "Guam";
case "gw":return "Guinea-Bissau";
case "gy":return "Guyana";

case "hk":return "Hong Kong";
case "hm":return "the Heard and McDonald Islands";
case "hn":return "Honduras";
case "hr":return "Croatia (Hrvatska)";
case "ht":return "Haiti";
case "hu":return "Hungary";

case "ic":return "the Canary Islands (reserved)";
case "id":return "Indonesia";
case "ie":return "Ireland";
case "il":return "Israel";
case "im":return "the Isle of Man";
case "in":return "India";
case "io":return "the British Indian Ocean Territory";
case "iq":return "Iraq";
case "ir":return "the Islamic Republic of Iran";
case "is":return "Iceland";
case "it":return "Italy";

case "je":return "Jersey";
case "jm":return "Jamaica";
case "jo":return "Jordan";
case "jp":return "Japan";

case "ke":return "Kenya";
case "kg":return "Kyrgyzstan";
case "kh":return "Cambodia";
case "ki":return "Kiribati";
case "km":return "Comoros";
case "kn":return "Saint Kitts and Nevis";
case "kp":return "the Democratic People's Republic of Korea";
case "kr":return "the Republic of Korea";
case "kw":return "Kuwait";
case "ky":return "the Cayman Islands";
case "kz":return "Kazakhstan";

case "la":return "the Lao People's Democratic Republic";
case "lb":return "Lebanon";
case "lc":return "Saint Lucia";
case "li":return "Liechtenstein";
case "lk":return "Sri Lanka";
case "lr":return "Liberia";
case "ls":return "Lesotho";
case "lt":return "Lithuania";
case "lu":return "Luxembourg";
case "lv":return "Latvia";
case "ly":return "the Libyan Arab Jamahiriya";

case "ma":return "Morocco";
case "mc":return "Monaco";
case "md":return "the Republic of Moldova";
case "me":return "Montenegro";
case "mf":return "Saint Martin (French part)";
case "mg":return "Madagascar";
case "mh":return "the Marshall Islands";
case "mk":return "the Former Yugoslav Republic of Macedonia";
case "ml":return "Mali";
case "mm":return "Myanmar";
case "mn":return "Mongolia";
case "mo":return "Macao";
case "mp":return "the Northern Mariana Islands";
case "mq":return "Martinique";
case "mr":return "Mauritania";
case "ms":return "Montserrat";
case "mt":return "Malta";
case "mu":return "Mauritius";
case "mv":return "the Maldives";
case "mw":return "Malawi";
case "mx":return "Mexico";
case "my":return "Malaysia";
case "mz":return "Mozambique";

case "na":return "Namibia";
case "nc":return "New Caledonia";
case "ne":return "Niger";
case "nf":return "Norfolk Island";
case "ng":return "Nigeria";
case "ni":return "Nicaragua";
case "nl":return "the Netherlands";
case "no":return "Norway";
case "np":return "Nepal";
case "nr":return "Nauru";
case "nt":return "the Saudi-Iraqi Neutral Zone (obsolete)";
case "nu":return "Niue";
case "nz":return "New Zealand";

case "om":return "Oman";

case "pa":return "Panama";
case "pe":return "Peru";
case "pf":return "French Polynesia";
case "pg":return "Papua New Guinea";
case "ph":return "the Philippines";
case "pk":return "Pakistan";
case "pl":return "Poland";
case "pm":return "Saint Pierre and Miquelon";
case "pn":return "Pitcairn";
case "pr":return "Puerto Rico";
case "ps":return "the Occupied Palestinian Territory";
case "pt":return "Portugal";
case "pw":return "Palau";
case "py":return "Paraguay";

case "qa":return "Qatar";
case "qo":return "Outlying Oceania";
case "qu":return "the European Union";

case "re":return "Re/union";
case "ro":return "Romania";
case "rs":return "Serbia";
case "ru":return "the Russian Federation";
case "rw":return "Rwanda";

case "sa":return "Saudi Arabia";
case "sb":return "the Solomon Islands";
case "sc":return "Seychelles";
case "sd":return "Sudan";
case "se":return "Sweden";
case "sf":return "Finland (obsolete)";
case "sg":return "Singapore";
case "sh":return "Saint Helena";
case "si":return "Slovenia";
case "sj":return "the Svalbard and Jan Mayen Islands";
case "sk":return "Slovakia (Slovak Republic)";
case "sl":return "Sierra Leone";
case "sm":return "San Marino";
case "sn":return "Senegal";
case "so":return "Somalia";
case "sr":return "Suriname";
case "st":return "Sa~o Tome/ and Pri/ncipe";
case "su":return "the Union of Soviet Socialist Republics (reserved)";
case "sv":return "El Salvador";
case "sy":return "the Syrian Arab Republic";
case "sz":return "Swaziland";

case "ta":return "Tristan da Cunha (reserved)";
case "tc":return "the Turks and Caicos Islands";
case "td":return "Chad";
case "tf":return "the French Southern Territories";
case "tg":return "Togo";
case "th":return "Thailand";
case "tj":return "Tajikistan";
case "tk":return "Tokelau";
case "tl":return "Timor-Leste";
case "tm":return "Turkmenistan";
case "tn":return "Tunisia";
case "to":return "Tonga";
case "tp":return "East Timor (obsolete)";
case "tr":return "Turkey";
case "tt":return "Trinidad and Tobago";
case "tv":return "Tuvalu";
case "tw":return "Taiwan Province of China";
case "tz":return "the United Republic of Tanzania";

case "ua":return "Ukraine";
case "ug":return "Uganda";
case "uk":return "the United Kingdom (reserved)";
case "um":return "the United States Minor Outlying Islands";
case "us":return "the United States of America";
case "uy":return "Uruguay";
case "uz":return "Uzbekistan";

case "va":return "the Holy See Vatican City State";
case "vc":return "Saint Vincent and the Grenadines";
case "ve":return "the Bolivarian Republic of Venezuela";
case "vg":return "the British Virgin Islands";
case "vi":return "the United States Virgin Islands";
case "vn":return "Viet Nam";
case "vu":return "Vanuatu";

case "wf":return "the Wallis and Futuna Islands";
case "ws":return "Samoa";

case "ye":return "Yemen";
case "yt":return "Mayotte";
case "yu":return "Yugoslavia (obsolete)";

case "xz":return "internatinal waters";

case "za":return "South Africa";
case "zm":return "Zambia";
case "zr":return "Zaire (obsolete)";
case "zw":return "Zimbabwe";
case "zz":return "unknown or invalid territory";
default:  return "Utopia";
}
}

/* info for *statename()
func: expand country & state names
args: string multi-letter ISO 3166-2 abbreviation
rtrn: string country & state name
*/
string *statename(string a){
string b,c,d,e,f;
a=lower_case(a);
if(sscanf(a,"%2s-%s",b,c)==2){
sscanf(a,"%2s-%s",b,c);
}else if(sscanf(a,"%2s",b)==1){
sscanf(a,"%2s",b);
}else{
return ({"a nameless country","somewhere","district"});
}
switch(b){
case "ad":d="Andorra";{
switch(c){
case "07":f="parish";e="Andorra la Vella";break;
case "02":f="parish";e="Canillo";break;
case "03":f="parish";e="Encamp";break;
case "08":f="parish";e="Escaldes-Engordany";break;
case "04":f="parish";e="La Massana";break;
case "05":f="parish";e="Ordino";break;
case "06":f="parish";e="Sant Julia\\ de Lo\\ria";break;
default:  f="parish";e="somewhere";break;
}
break;}
case "au":d="Australia";{
switch(c){
case "nsw":f="state";e="New South Wales";break;
case "qld":f="state";e="Queensland";break;
case "sa": f="state";e="South Australia";break;
case "tas":f="state";e="Tasmania";break;
case "vic":f="state";e="Victoria";break;
case "wa": f="state";e="Western Australia";break;
case "act":f="territory";e="the Australian Capital Territory";break;
case "nt": f="territory";e="the Northern Territory";break;
case "cc": f="external territory";e="the Cocos (Keeling) Islands";break;
case "cx": f="external territory";e="Christmas Island";break;
case "hm": f="external territory";e="Heard Island and McDonald Islands";break;
case "nf": f="external territory";e="Norfolk Island";break;
default:   f="district";e="somewhere";break;
}
break;}
case "ca":d="Canada";{
switch(c){
case "ab":f="province";e="Alberta";break;
case "bc":f="province";e="British Columbia";break;
case "mb":f="province";e="Manitoba";break;
case "nb":f="province";e="New Brunswick";break;
case "nl":f="province";e="Newfoundland and Labrador";break;
case "ns":f="province";e="Nova Scotia";break;
case "on":f="province";e="Ontario";break;
case "pe":f="province";e="Prince Edward Island";break;
case "qc":f="province";e="Quebec";break;
case "sk":f="province";e="Saskatchewan";break;
case "nt":f="territory";e="the Northwest Territories";break;
case "nu":f="territory";e="Nunavut";break;
case "yt":f="territory";e="the Yukon";break;
default:  f="district";e="somewhere";break;
}
break;}
case "de":d="Germany";{
switch(d){
case "bw":f="land";e="Baden-Wu\"ttemberg";break;
case "by":f="land";e="Bayern";break;
case "be":f="land";e="Berlin";break;
case "bb":f="land";e="Brandenburg";break;
case "hb":f="land";e="Bremen";break;
case "hh":f="land";e="Hamburg";break;
case "he":f="land";e="Hessen";break;
case "mv":f="land";e="Mecklenburg-Vorpommern";break;
case "ni":f="land";e="Niedersachsen";break;
case "nw":f="land";e="Nordrhein-Westfalen";break;
case "rp":f="land";e="Rheinland-Pfalz";break;
case "sl":f="land";e="Saarland";break;
case "sn":f="land";e="Sachsen";break;
case "st":f="land";e="Sachsen-Anhalt";break;
case "sh":f="land";e="Schleswig-Holstein";break;
case "th":f="land";e="Thu\"ringen";break;
default:  f="district";e="somewhere";break;
}
break;}
case "us":d="the United States of America";{
switch(c){
case "al":f="state";e="Alabama";break;
case "ak":f="state";e="Alaska";break;
case "az":f="state";e="Arizona";break;
case "ar":f="state";e="Arkansas";break;
case "ca":f="state";e="California";break;
case "co":f="state";e="Colorado";break;
case "ct":f="state";e="Connecticut";break;
case "de":f="state";e="Delaware";break;
case "fl":f="state";e="Florida";break;
case "ga":f="state";e="Georgia";break;
case "hi":f="state";e="Hawaii";break;
case "id":f="state";e="Idaho";break;
case "il":f="state";e="Illinois";break;
case "in":f="state";e="Indiana";break;
case "ia":f="state";e="Iowa";break;
case "ks":f="state";e="Kansas";break;
case "ky":f="state";e="Kentucky";break;
case "la":f="state";e="Louisiana";break;
case "me":f="state";e="Maine";break;
case "md":f="state";e="Maryland";break;
case "ma":f="state";e="Massachussets";break;
case "mi":f="state";e="Michigan";break;
case "mn":f="state";e="Minnesota";break;
case "ms":f="state";e="Mississippi";break;
case "mo":f="state";e="Missouri";break;
case "mt":f="state";e="Montana";break;
case "ne":f="state";e="Nebraska";break;
case "nv":f="state";e="Nevada";break;
case "nh":f="state";e="New Hampshire";break;
case "nj":f="state";e="New Jersey";break;
case "nm":f="state";e="New Mexico";break;
case "ny":f="state";e="New York";break;
case "nc":f="state";e="North Carolina";break;
case "nd":f="state";e="North Dakota";break;
case "oh":f="state";e="Ohio";break;
case "ok":f="state";e="Oklahoma";break;
case "or":f="state";e="Oregon";break;
case "pa":f="state";e="Pennsylvania";break;
case "ri":f="state";e="Rhode Island and Providence Plantations";break;
case "sc":f="state";e="South Carolina";break;
case "sd":f="state";e="South Dakota";break;
case "tn":f="state";e="Tennessee";break;
case "tx":f="state";e="Texas";break;
case "ut":f="state";e="Utah";break;
case "vt":f="state";e="Vermont";break;
case "va":f="state";e="Virginia";break;
case "wa":f="state";e="Washington";break;
case "wv":f="state";e="West Virginia";break;
case "wi":f="state";e="Wisconsin";break;
case "wy":f="state";e="Wyoming";break;
case "dc":f="District";e="Columbia";break;
case "as":f="outlying territory";e="American Samoa";break;
case "gu":f="outlying territory";e="Guam";break;
case "mp":f="outlying territory";e="the Northern Mariana Islands";break;
case "pr":f="outlying territory";e="Puerto Rico";break;
case "vi":f="outlying territory";e="the US Virgin Islands";break;
case "um":f="outlying territory (obsolete)";e="the US minor outlying islands";break;
case "fm":f="outlying territory (obsolete)";e="the Federated States of Micronesia";break;
case "mh":f="outlying territory (obsolete)";e="the Marshall Islands";break;
case "pw":f="outlying territory (obsolete)";e="Palau";break;
default:  f="district";e="somewhere";break;
}
break;}
default:d="a nameless country";e="somewhere";f="district";break;
}
return ({d,e,f});
}

/* info for city_loci()
func: get locations from city names
args: string iso 3166-2 state, string county, string city name
rtrn: float array ({latitude,longitude,altitude (m)})
*/
float *city_loci(string state,string county,string city){
switch(lower_case(state)){
case "us-ia":
switch(lower_case(county)){
case "henry county":
switch(lower_case(city)){
case "coppock":         return ({41.1626,-091.7150,  193.0});
case "hillsboro":       return ({40.8363,-091.7136,  223.0});
case "mount pleasant":  return ({40.9670,-091.5512,  218.0});
case "mount union":     return ({41.0575,-091.3908,  222.0});
case "new london":      return ({40.9252,-091.4011,  232.0});
case "olds":            return ({41.1333,-091.5456,  222.0});
case "rome":            return ({40.9830,-091.6830,  182.0});
case "salem":           return ({40.8531,-091.6208,  218.0});
case "swedesburg":      return ({41.1053,-091.5472,  223.0});
case "wayland":         return ({41.1476,-091.6597,  223.0});
case "westwood":        return ({40.9642,-091.6283,  213.0});
case "winfield":        return ({41.1271,-091.4374,  210.0});
default:                return ({41.0000,-091.5000,  200.0});
}
default:                     return({41.5000,-093.5000,  335.0});
}
case "us-fl":
switch(lower_case(county)){
case "brevard county":
switch(lower_case(city)){
case "merritt island":   return ({28.3578,-080.6847,   1.0});
case "rockledge":        return ({28.3250,-080.7328,   7.0});
default:                 return ({00.0000, 000.0000,   0.0});
}
case "orange county":
switch(lower_case(city)){
case "apopka":           return ({28.6761,-081.5106,  40.0});
case "belle isle":       return ({28.4644,-081.3503,  29.0});
case "eatonville":       return ({28.6186,-081.3833,  29.0});
case "edgewood":         return ({28.4869,-081.3733,  29.0});
case "maitland":         return ({28.6269,-081.3669,  27.0});
case "oakland":          return ({28.5544,-081.6311,  37.0});
case "ocoee":            return ({28.5742,-081.5306,  37.0});
case "orlando":          return ({28.5436,-081.3728,  34.0});
case "plymouth":         return ({28.6761,-081.5106,  40.0});
case "windmere":         return ({28.4967,-081.5339,  37.0});
case "winter garden":    return ({28.5603,-081.5842,  38.0});
case "winter park":      return ({28.5961,-081.3467,  28.0});
case "bay lake":         return ({28.3914,-081.5667,  28.0});
case "lake buena vista": return ({28.3775,-081.5217,  29.0});
default:                 return ({28.5100,-081.3200,   0.0});
}
case "seminole county":
switch(lower_case(city)){
case "altamonte springs":return ({28.6614,-081.3919,  26.0});
case "casselberry":      return ({28.6603,-081.3233,  17.0});
case "lake mary":        return ({28.7575,-081.3292,  18.9});
case "longwood":         return ({28.7019,-081.3450,  22.9});
case "oviedo":           return ({28.6597,-081.1958,  14.6});
case "sanford":          return ({28.8003,-081.2733,  11.0});
case "winter springs":   return ({28.6933,-081.2789,  16.0});
case "forest city":      return ({28.6628,-081.4428,  29.0});
case "heathrow":         return ({28.7728,-081.3714,  15.0});
case "tuskawilla":       return ({28.6933,-081.2789,  16.0});
case "chuluota":         return ({28.6414,-081.1244,  18.0});
case "fern park":        return ({28.6481,-081.3464,  29.0});
case "geneva":           return ({28.7372,-081.1178,  19.0});
case "midway":           return ({28.7897,-081.2289,  25.0});
case "wekiwa springs":   return ({28.6969,-081.4258,  16.0});
case "winter park":      return ({28.5961,-081.3467,  28.0});
default:                 return ({28.7100,-081.2300,   0.0});
}
default:                     return ({28.1333,-081.6317,   0.0});
}
case "us-il":
switch(lower_case(county)){
case "kane county":
switch(lower_case(city)){
case "aurora":           return ({41.7633,-088.2900, 220.0});
case "batavia":          return ({41.8489,-088.3083, 203.0});
case "big rock":         return ({41.7639,-088.5469, 216.0});
case "burlington":       return ({42.0522,-088.5483, 280.0});
case "campton hills":    return ({41.9339,-088.4011, 271.0});
case "east dundee":      return ({42.1018,-088.2702, 249.0});
case "elgin":            return ({42.0394,-088.2886, 250.0});
case "geneva":           return ({41.8761,-088.3219, 224.0});
case "hampshire":        return ({42.0989,-088.5258, 269.0});
case "kaneville":        return ({41.8353,-088.5219, 259.0});
case "st. charles":      return ({41.9137,-088.3109, 222.0});
case "sugar grove":      return ({41.7725,-088.4425, 213.0});
case "virgil":           return ({41.9567,-088.5286, 264.0});
default:                 return ({41.9500,-088.4300, 200.0});
}
default:                     return ({40.0000,-089.0000,   0.0});
}
case "us-nm":
switch(lower_case(county)){
case "bernalillo county":
switch(lower_case(city)){
case "albuquerque":      return ({35.1108,-106.6100,1619.1});
case "alburquerque":     return ({35.1108,-106.6100,1619.1});
default:                 return ({35.0500,-106.6700,   0.0});
}
default:                     return ({34.0000, 106.0000,   0.0});
}
default:                         return ({ 0.0000,   0.0000,   0.0});
}
}
/* EOF */

Code: [Select]
/* ** HEADER for ../funs/astro/stars.c **
2005-Jul-17 - T. Cook
star data
-----
2011-Jul-17 - T. Cook began altering for DeadSouls
2010-Jun-01 - T. Cook reformatted for notepad++
2010-Jan-03 - T. Cook switched back to a function
2009-Mar-15 - T. Cook simplified to single #declare, added names, absolute magnitude, distance
2005-Sep-28 - T. Cook tidied up
to do:
add more stars (limit 1000?)
*/
mixed *STARS(){
return ({
// RA     dec    v mag  a mag  dist   abbr     name
({ 06.45, -16.7, -1.46,  1.43,    9, "CMa A", "Sirius" }),
({ 06.24, -52.7, -0.73, -5.64,  310, "Car A", "Canopus" }),
({ 14.40, -60.8, -0.29,  4.06,    4, "Cen A", "Rigil Kentaurus" }),
({ 14.16,  19.2, -0.05, -0.31,   37, "Boo A", "Arcturus" }),
({ 18.37,  38.8,  0.03,  0.58,   25, "Lyr A", "Vega" }),
({ 05.17,  46.0,  0.07, -0.49,   42, "Aur A", "Capella" }),
({ 05.15, -08.2,  0.15, -6.72,  770, "Ori B", "Rigel" }),
({ 07.39,  05.2,  0.36,  2.64,   11, "CMi A", "Procyon" }),
({ 01.38, -57.2,  0.45, -2.77,  144, "Eri A", "Achernar" }),
({ 05.55,  07.4,  0.55, -5.04,  430, "Ori A", "Betelgeuse" }),
({ 14.04, -60.4,  0.61, -5.42,  530, "Cen B", "Hadar" }),
({ 19.51,  08.9,  0.77,  2.21,   17, "Aql A", "Altair" }),
({ 12.27, -63.1,  0.79, -4.17,  320, "Crx A", "Acrux" }),
({ 04.36,  16.5,  0.86, -0.64,   65, "Tau A", "Aldebaran" }),
({ 16.29, -26.4,  0.95, -5.39,  600, "Sco A", "Antares" }),
({ 13.25, -11.2,  0.97, -3.56,  260, "Vir A", "Spica" }),
({ 07.45,  28.0,  1.14,  1.07,   34, "Gem B", "Pollux" }),
({ 22.58, -29.6,  1.15,  1.72,   25, "PsA A", "Fomalhaut" }),
({ 20.41,  45.3,  1.24, -8.75, 3000, "Cyg A", "Deneb" }),
({ 12.48, -59.7,  1.26, -3.91,  350, "Crx B", "Mimosa" }),
({ 10.08,  12.0,  1.36, -0.52,   78, "Leo A", "Regulus" }),
({ 06.59, -29.0,  1.50, -4.10,  430, "CMa E", "Adhara" }),
({ 07.35,  31.9,  1.58,  0.59,   52, "Gem A", "Castor" }),
({ 17.34, -37.1,  1.62, -5.05,  700, "Sco L", "Shaula" }),
({ 12.31, -57.1,  1.63, -0.52,   88, "Crx G", "Gacrux" }),
({ 05.25,  06.3,  1.64, -2.72,  240, "Ori G", "Bellatrix" }),
({ 05.26,  28.6,  1.66, -1.36,  130, "Tau B", "Elnath" }),
({ 09.13, -69.7,  1.67, -0.99,  111, "Car B", "Miaplacidus" }),
({ 05.36, -01.2,  1.69, -6.38, 1300, "Ori E", "Alnilam" }),
({ 22.08, -47.0,  1.74, -0.72,  101, "Gru A", "Alnair" }),
({ 05.41, -01.9,  1.75, -5.25,  820, "Ori Z", "Alnitak" }),
({ 12.54,  56.0,  1.77, -0.20,   81, "UMa E", "Alioth" }),
({ 03.24,  49.9,  1.80, -4.49,  590, "Per A", "Mirfak" }),
({ 11.04,  61.8,  1.80, -1.09,  124, "UMa A", "Dubhe" }),
({ 08.10, -47.3,  1.81, -5.25,  840, "Vel G", "Regor" }),
({ 07.08, -26.4,  1.83, -6.87, 1800, "CMa D", "Wezen" }),
({ 18.24, -34.4,  1.84, -1.39,  145, "Sag E", "Kaus Australis" }),
({ 13.48,  49.3,  1.86, -0.59,  101, "UMa H", "Alkaid" }),
({ 17.37, -43.0,  1.86, -2.75,  270, "Sco t", "Sargas" }),
({ 08.23, -59.5,  1.87, -4.57,  630, "Car E", "Avior" }),
({ 06.00,  44.9,  1.90, -0.10,   82, "Aur B", "Menkalinan" }),
({ 16.49, -69.0,  1.92, -3.61,  420, "TrA A", "Atria" }),
({ 06.38,  16.4,  1.93, -0.60,  105, "Gem G", "Alhena" }),
({ 20.26, -56.7,  1.93, -1.82,  180, "Pav A", "Peacock" }),
({ 08.45, -54.7,  1.95,  0.01,   80, "Vel D", "Koo She" }),
({ 06.23, -18.0,  1.98, -3.95,  500, "CMa B", "Mirzam" }),
({ 09.28, -08.7,  1.98, -1.70,  180, "Hyd A", "Alphard" }),
({ 02.32,  89.3,  1.99, -3.62,  430, "UMi A", "Polaris" }),
({ 10.20,  19.8,  2.00, -0.93,  126, "Leo G", "Algieba" }),
({ 02.07,  23.5,  2.01,  0.48,   66, "Ari A", "Hamal" }),
({ 00.44, -18.0,  2.04, -0.30,   96, "Cet B", "Diphda" }),
({ 18.55, -26.3,  2.05, -2.14,  220, "Sag S", "Nunki" }),
({ 14.07, -36.4,  2.06,  0.70,   61, "Cen t", "Menkent" }),
({ 00.08,  29.1,  2.07, -0.30,   97, "And A", "Alpheratz" }),
({ 01.10,  35.6,  2.07, -1.86,  200, "And B", "Mirach" }),
({ 05.48, -09.7,  2.07, -4.65,  720, "Ori K", "Saiph" }),
({ 14.51,  74.2,  2.07, -0.87,  127, "UMi B", "Kochab" }),
({ 22.43, -46.9,  2.07, -1.52,  170, "Gru B", "Al Dhanab" }),
({ 17.35,  12.6,  2.08,  1.30,   47, "Oph A", "Rasalhague" }),
({ 03.08,  41.0,  2.09, -0.18,   93, "Per B", "Algol" }),
({ 02.04,  42.3,  2.10, -3.08,  360, "And G", "Almach" }),
({ 11.49,  14.6,  2.14,  1.92,   36, "Leo B", "Denebola" }),
({ 00.57,  60.7,  2.15, -4.22,  610, "Cas G", "Cih" }),
({ 12.42, -49.0,  2.20, -0.81,  130, "Cen G", "Muhlifain" }),
({ 08.04, -40.0,  2.21, -5.95, 1400, "Pup Z", "Naos" }),
({ 09.17, -59.3,  2.21, -4.42,  690, "Car I", "Aspidiske" }),
({ 15.35,  26.7,  2.22,  0.42,   75, "CrB A", "Alphecca" }),
({ 09.08, -43.4,  2.23, -3.99,  570, "Vel L", "Suhail" }),
({ 13.24,  54.9,  2.23,  0.33,   78, "UMa Z", "Mizar" }),
({ 20.22,  40.3,  2.23, -6.12, 1500, "Cyg G", "Sadr" }),
({ 00.41,  56.5,  2.24, -1.99,  230, "Cas A", "Schedar" }),
({ 17.57,  51.5,  2.24, -1.04,  148, "Dra G", "Eltanin" }),
({ 05.32, -00.3,  2.25, -4.99,  920, "Ori D", "Mintaka" }),
({ 00.09,  59.2,  2.28,  1.17,   55, "Cas B", "Caph" }),
({ 13.40, -53.5,  2.29, -3.02,  380, "Cen E", "" }),
({ 16.00, -22.6,  2.29, -3.16,  400, "Sco D", "Dschubba" }),
({ 16.50, -34.3,  2.29,  0.78,   65, "Sco E", "Wei" }),
({ 14.42, -47.4,  2.30, -3.83,  550, "Lup A", "Men" }),
({ 14.36, -42.2,  2.33, -2.55,  310, "Cen E", "" }),
({ 11.02,  56.4,  2.34,  0.41,   79, "UMa B", "Merak" }),
({ 14.45,  27.1,  2.35, -1.69,  210, "Boo E", "Izar" }),
({ 21.44,  09.9,  2.38, -4.19,  670, "Peg E", "Enif" }),
({ 17.42, -39.0,  2.39, -3.38,  460, "Sco K", "Girtab" }),
({ 00.26, -42.3,  2.40,  0.52,   77, "Phe A", "Ankaa" }),
({ 11.54,  53.7,  2.41,  0.36,   84, "UMa G", "Phecda" }),
({ 17.10, -15.7,  2.43,  0.37,   84, "Oph E", "Sabik" }),
({ 23.04,  28.1,  2.44, -1.49,  200, "Peg B", "Scheat" }),
({ 07.24, -29.3,  2.45, -7.51, 3000, "CMa E", "Aludra" }),
({ 21.19,  62.6,  2.45,  1.58,   49, "Cep A", "Alderamin" }),
({ 09.22, -55.0,  2.47, -3.62,  540, "Vel K", "Markeb" }),
({ 20.46,  34.0,  2.48,  0.76,   72, "Cyg E", "Gienah" }),
({ 23.05,  15.2,  2.49, -0.67,  140, "Peg A", "Markab" }),
({ 03.02,  04.1,  2.54, -1.61,  220, "Cet A", "Menkar" }),
({ 16.37, -10.6,  2.54, -3.20,  460, "Oph Z", "Han" }),
({ 13.56, -47.3,  2.55, -2.81,  390, "Cen Z", "Al Nair al Kentaurus" }),
({ 11.14,  20.5,  2.56,  1.32,   58, "Leo D", "Zosma" }),
({ 16.05, -19.8,  2.56, -3.50,  530, "Sco B", "Graffias" }),
({ 05.33, -17.8,  2.58, -5.40, 1300, "Lep A", "Arneb" }),
({ 12.08, -50.7,  2.58, -2.84,  400, "Cen D", "" }),
({ 12.16, -17.5,  2.58, -0.94,  165, "Cor G", "Gienah Ghurab" })
});
}
/* EOF */

Offline z993126

  • BFF
  • ***
  • Posts: 128
    • View Profile
Re: Sky map and astronomical functions
« Reply #3 on: July 23, 2011, 01:25:56 AM »
Code: [Select]
/* ** HEADER for ../funs/astro/earth_sun.c **
2005-May-27 - T. Cook
functions returning various solar position elements
sun position elements from aa.usno.navy.mil
sun position from users.zoominternet.net/~matto
-----
2011-Jul-21 - T. Cook removed unneeded functions
2011-Jul-17 - T. Cook began altering for DeadSouls
2011-Jul-16 - T. Cook added calcSunRSJD() and calcSolNoonJD()
2011-Jul-12 - T. Cook tidied up
2010-Jul-06 - T. Cook added calc[Last/Next][Sunrise/Sunset]()
2010-Jul-04 - T. Cook corrected calcSunRS* issues
2010-Jun-11 - T. Cook fixed RiseSetNew() to use N=000, E=090
2010-Jun-01 - T. Cook reformatted for notepad++
2010-Jan-24 - T. Cook added matto's sunrise/set function
2009-Dec-29 - T. Cook modified for tublib
2005-Sep-28 - T. Cook tidied up
to do:
simplify/streamline things! (combine functions?)
remove unneeded UTC & Local functions
*/
#pragma strong_types

//------------------------------------------------------------------//
#include "/realms/ardneh/area/customdefs.h"
#include "../../inc/astro.h"
#include "../../inc/math.h"

inherit MY_FUNS "/math.c";
inherit MY_FUNS "/astro/chronfuns.c";
inherit MY_FUNS "/lang/longnumbers.c";
inherit MY_FUNS "/geo/timezone.c";
//------------------------------------------------------------------//
float calcGeomMeanLongSun( float f_t );
float calcGeomMeanAnomalySun( float f_t );
float calcEccentricityEarthOrbit( float f_t );
float calcSunEqOfCenter( float f_t );
float calcSunTrueLong( float f_t );
float calcSunTrueAnomaly( float f_t );
float calcSunRadVector( float f_t );
float calcSunApparentLong( float f_t );
float calcMeanObliquityOfEcliptic( float f_t );
float calcObliquityCorrection( float f_t );
float calcSunRtAscension( float f_t );
float calcSunDeclination( float f_t );
float calcEquationOfTime( float f_t );
float calcHourAngleSunRS( float f_latitude, float f_solarDec, int i_RS );
mixed *calcSunRSJD( int i_JD, float f_dayfrac, float f_latitude, float f_longitude, string s_type, int i_RS );
mixed *calcSolNoonJD( int i_JD, float f_longitude );
float *calcSunPos( int i_JD, float f_df);
/********************************************************************/
/* info for calcGeomMeanLongSun()
func: calculate Geometric Mean Longitude of the Sun (degrees)
args: float Julian centuries since J2000.0
rtrn: float Geometric Mean Longitude of the Sun in degrees
*/
float calcGeomMeanLongSun( float f_t ){
return fwrap( 280.46646 + f_t * ( 36000.76983 + 0.0003032 * f_t ), 0.0, 359.99999 );
}
/* info for calcGeomMeanAnomalySun()
func: calculate Geometric Mean Anomaly of the Sun (degrees)
args: Julian centuries since J2000.0
rtrn: Geometric Mean Anomaly of the Sun in degrees
*/
float calcGeomMeanAnomalySun( float f_t ){
return 357.52911 + f_t * ( 35999.05029 - 0.0001537 * f_t );
}
/* info for calcEccentricityEarthOrbit()
func: calculate the eccentricity of Earth's orbit (unitless)
args: Julian centuries since J2000.0
rtrn: eccentricity of Earth's orbit
*/
float calcEccentricityEarthOrbit( float f_t ){
return 0.016708634 - f_t * ( 0.000042037 + 0.0000001267 * f_t );
}
/* info for calcSunEqOfCenter()
func: calculate the equation of centre for the Sun (degrees)
args: Julian centuries since J2000.0
rtrn: equation of centre for the Sun in degrees
*/
float calcSunEqOfCenter( float f_t ){
float f_mrad = deg2rad( calcGeomMeanAnomalySun( f_t ) ),
f_sinm = sin( f_mrad ),
f_sin2m = sin( 2 * f_mrad ),
f_sin3m = sin( 3 * f_mrad );
return f_sinm * ( 1.914602 - f_t * ( 0.004817 + 0.000014 * f_t ) ) + f_sin2m * ( 0.019993 - 0.000101 * f_t ) + f_sin3m * 0.000289;
}
/* info for calcSunTrueLong()
func: calculate the true longitude of the Sun (degrees)
args: Julian centuries since J2000.0
rtrn: true longitude of the Sun in degrees
*/
float calcSunTrueLong( float f_t ){
return calcGeomMeanLongSun( f_t) + calcSunEqOfCenter( f_t );
}
/* info for calcSunTrueAnomaly()
func: calculate the true anomaly of the Sun (degrees)
args: Julian centuries since J2000.0
rtrn: true anomaly of the Sun in degrees
*/
float calcSunTrueAnomaly( float f_t ){
return calcGeomMeanAnomalySun( f_t ) + calcSunEqOfCenter( f_t );
}
/* info for calcSunRadVector()
func: calculate the Earth's distance from the Sun in (AUs)
args: Julian centuries since J2000.0
rtrn: Sun radius vector in AUs
*/
float calcSunRadVector( float f_t ){
float f_e = calcEccentricityEarthOrbit( f_t );
return ( 1.000001018 * ( 1 - f_e * f_e ) ) / ( 1 + f_e * cos( deg2rad( calcSunTrueAnomaly( f_t ) ) ) );
}
/* info for calcSunApparentLong()
func: calculate the apparent longitude of the Sun (degrees)
args: Julian centuries since J2000.0
rtrn: Sun's apparent longitude in degrees
*/
float calcSunApparentLong( float f_t ){
return calcSunTrueLong( f_t ) - 0.00569 - 0.00478 * sin( deg2rad( 125.04 - 1934.136 * f_t ) );
}
/* info for calcMeanObliquityOfEcliptic()
func: calculate the mean obliquity of the ecliptic (degrees)
args: Julian centuries since J2000.0
rtrn: mean obliquity of ecliptic in degrees
*/
float calcMeanObliquityOfEcliptic( float f_t ){
return 23.0 + ( 26.0 + ( ( 21.448 - f_t * ( 46.8150 + f_t * ( 0.00059 - f_t * ( 0.001813 ) ) ) ) / 60.0 ) ) / 60.0;
}
/* info for calcObliquityCorrection()
func: calculate the corrected obliquity of the ecliptic (degrees)
args: Julian centuries since J2000.0
rtrn: corrected obliquity of ecliptic in degrees
*/
float calcObliquityCorrection( float f_t ){
return calcMeanObliquityOfEcliptic( f_t ) + 0.00256 * cos( deg2rad( 125.04 - 1934.136 * f_t ) );
}
/* info for calcSunRtAscension()
func: calculate the right ascension of the Sun (degrees)
args: Julian centuries since J2000.0
rtrn: Sun's right ascension in degrees
*/
float calcSunRtAscension( float f_t ){
float f_lambda = calcSunApparentLong( f_t );
return rad2deg( atanb( ( cos( deg2rad( calcObliquityCorrection( f_t ) ) ) *
sin( deg2rad( f_lambda ) ) ), ( cos( deg2rad( f_lambda ) ) ) ) );
}
/* info for calcSunDeclination()
func: calculate the declination of the Sun (degrees)
args: Julian centuries since J2000.0
rtrn: Sun's declination in degrees
*/
float calcSunDeclination( float f_t ){
//float e=calcObliquityCorrection(t);
//float lambda=calcSunApparentLong(t);
//float sint=sin(deg2rad(e)) * sin(deg2rad(lambda));
//float theta=rad2deg(asin(sint));
//return theta;
return rad2deg( asin( sin( deg2rad( calcObliquityCorrection( f_t ) ) ) * sin( deg2rad( calcSunApparentLong( f_t ) ) ) ) );
}
/* info for calcEquationOfTime()
func: calculate the difference between true and mean solar times
args: Julian centuries since J2000.0
rtrn: equation of time in minutes of time
*/
float calcEquationOfTime( float f_t ){
float f_epsilon, f_l0, f_e, f_m, f_y, f_sin2l0, f_sinm, f_cos2l0, f_sin4l0, f_sin2m, f_Etime;
f_epsilon=calcObliquityCorrection( f_t );
f_l0 = calcGeomMeanLongSun( f_t );
f_e = calcEccentricityEarthOrbit( f_t );
f_m = calcGeomMeanAnomalySun( f_t );
f_y = tan( deg2rad( f_epsilon ) / 2.0 );
f_y *= f_y;
f_sin2l0 = sin( 2.0 * deg2rad( f_l0 ) );
f_sinm = sin( deg2rad( f_m ) );
f_cos2l0 = cos( 2.0 * deg2rad( f_l0 ) );
f_sin4l0 = sin( 4.0 * deg2rad( f_l0 ) );
f_sin2m = sin( 2.0 * deg2rad( f_m ) );
f_Etime = f_y * f_sin2l0 - 2.0 * f_e * f_sinm + 4.0 * f_e * f_y * f_sinm *
f_cos2l0 - 0.5 * f_y * f_y * f_sin4l0 - 1.25 * f_e * f_e * f_sin2m;
return rad2deg( f_Etime ) * 4.0;
}
/* info for calcHourAngleRS()
func: calculate the sunrise hour angle (radians)
args: float latitude in degrees, float solar declination in degrees, int +1 for rise -1 for set
rtrn: float hour angle of sunrise in radians
*/
float calcHourAngleSunRS( float f_latitude, float f_solarDec, int i_RS ){
float f_latRad, f_sdRad, f_HAarg, f_HA;
if( f_latitude > 89.8 ){ f_latitude = 89.8; }
else if( f_latitude < -89.8 ){ f_latitude = -89.8; }
f_latRad = deg2rad( f_latitude );
f_sdRad = deg2rad( f_solarDec );
f_HAarg = ( cos( deg2rad( 90.833 ) ) / ( cos( f_latRad ) * cos( f_sdRad ) ) - tan( f_latRad ) * tan( f_sdRad ) );
if( catch( acos( cos( deg2rad( 90.833 ) ) / ( cos( f_latRad) * cos( f_sdRad ) ) - tan( f_latRad ) * tan( f_sdRad ) ) ) ){
return -1.0; // acos returns values 0..PI
}else{
f_HA=( acos( cos( deg2rad( 90.833 ) ) / ( cos( f_latRad ) * cos( f_sdRad ) ) - tan( f_latRad ) * tan( f_sdRad ) ) );
}
return f_HA * i_RS;
}

/* info for calcSunRSJD()
func: calculate the Julian Date of sunrise, set, or twilight for given day at the given location on Earth
args: Julian day, latitude, longitude, type, RS (rise=1, set=-1)
rtrn: mixed array sunrise/set/twilight time ({ int Julian date, float dayfrac })
*/
mixed *calcSunRSJD( int i_JD, float f_dayfrac, float f_latitude, float f_longitude, string s_type, int i_RS ){
mixed *ma_tnoon;
float f_t, f_tnoon, f_eqTime, f_solarDec, f_hourAngle;
float f_delta, f_timeDiff, f_timeUTC, f_newt;
string s_season;
if( f_latitude > 89.8 ){ f_latitude = 89.8; }
else if( f_latitude < -89.8 ){ f_latitude = -89.8; }
ma_tnoon = calcSolNoonJD( i_JD, f_longitude );
f_t = calcTimeJulianCent( i_JD + f_dayfrac );
f_tnoon = calcTimeJulianCent( ma_tnoon[0] + ma_tnoon[1] );
f_eqTime = calcEquationOfTime( f_tnoon );
f_solarDec = calcSunDeclination( f_tnoon );
f_hourAngle = calcHourAngleSunRS( f_latitude, f_solarDec, -i_RS);
f_delta = f_longitude - rad2deg( f_hourAngle );
f_timeDiff = 4 * f_delta;
f_timeUTC = 720 - f_timeDiff - f_eqTime;
f_newt = calcTimeJulianCent( i_JD + f_timeUTC / 1440.0 );
f_eqTime = calcEquationOfTime( f_newt );
f_solarDec = calcSunDeclination( f_newt );
f_hourAngle = calcHourAngleSunRS( f_latitude, f_solarDec, -i_RS );
// catch polar region rise/set information; find last or next event
s_season = season( i_JD, f_dayfrac, f_latitude )[0];
while( f_hourAngle == -1.0 ){
if(
( f_latitude >  66.4 && ( s_season == "spring" || s_season == "summer" || s_season == "midnight sun" ) ) ||
( f_latitude < -66.4 && ( s_season == "autumn" || s_season == "winter" || s_season == "midnight sun" ) )
){
if( i_RS == 1 ){ --i_JD; }else if( i_RS == -1 ){ ++i_JD; }else{ return ({ -1, 0.0 }); }
}else if(
( f_latitude < -66.4 && ( s_season == "spring" || s_season == "summer" || s_season == "polar night" ) ) ||
( f_latitude >  66.4 && ( s_season == "autumn" || s_season == "winter" || s_season == "polar night" ) )
){
if( i_RS == 1 ){ ++i_JD; }else if( i_RS == -1 ){ --i_JD; }else{ return ({ -1, 0.0 }); }
}else{
return ({ -1, 0.0 });
}
f_delta = f_longitude - rad2deg( f_hourAngle );
f_timeDiff = 4 * f_delta;
f_timeUTC = 720 - f_timeDiff - f_eqTime;
f_newt = calcTimeJulianCent( i_JD + f_timeUTC / 1440.0 );
f_eqTime = calcEquationOfTime( f_newt );
f_solarDec = calcSunDeclination( f_newt );
f_hourAngle = calcHourAngleSunRS( f_latitude, f_solarDec, -i_RS );
}
switch( s_type ){
case "civil":   f_delta = f_longitude - rad2deg( f_hourAngle ) +  6 * i_RS; break;
case "nautical":f_delta = f_longitude - rad2deg( f_hourAngle ) + 12 * i_RS; break;
case "astro":   f_delta = f_longitude - rad2deg( f_hourAngle ) + 18 * i_RS; break;
default:        f_delta = f_longitude - rad2deg( f_hourAngle );             break;
}
f_timeDiff = 4 * f_delta;
f_dayfrac = 720 - f_timeDiff - f_eqTime;
f_dayfrac = f_dayfrac / 1440.0 - 0.5;
while( f_dayfrac > 1.0 ){ ++i_JD; --f_dayfrac; }
while( f_dayfrac < 0.0 ){ --i_JD; ++f_dayfrac; }
return ({ i_JD, f_dayfrac });
}
/* info for *calcSolNoonJD()
func: calculate the Julian Date of solar noon for given day at given longitude on Earth
args: int Julian Date and float longitude
rtrn: mixed array noon time in ({ Julian Date, day fraction })
*/
mixed *calcSolNoonJD( int i_JD, float f_longitude ){
float f_JD = to_float( i_JD );
float f_dayfrac = 720 - ( f_longitude * 4 ) - calcEquationOfTime( calcTimeJulianCent( f_JD + 0.5 + f_longitude / 360.0 ) );
f_dayfrac = f_dayfrac / 1440.0 - 0.5;
if( f_dayfrac > 1.0 ){ ++i_JD; --f_dayfrac; }
if( f_dayfrac < 0.0 ){ --i_JD; ++f_dayfrac; }
return ({ i_JD, f_dayfrac });
}

/* info for *calcSunPos()
func: calculate apparent position of Sol from centre of Earth
args: Julian day, dayfrac
rtrn: array ({ RA in decimal hours, decl in decimal degrees, longitude in degrees, obliquity, distance in AU })
*/
float *calcSunPos( int i_JD, float f_dayfrac ){
// variable declarations
float f_J2, f_ET, f_ET2, f_VP, f_P, f_M0, f_MN, f_T0, f_S, f_AM, f_V, f_R, f_L, f_Z, f_OB, f_DC, f_AR;
f_J2 = i_JD + f_dayfrac;
f_ET  = 0.016718;
f_ET2 = f_ET * f_ET;
f_VP  = 0.0000822;
f_P   = 4.93204;
f_M0  = 2.12344;
f_MN  = 0.0172019;
f_T0  = 2444000.5;
f_S   = 2415020.5;
f_P += ( f_J2 - f_T0 ) * f_VP / 100;
f_AM = f_M0 + f_MN * ( f_J2 - f_T0 );
f_AM -= 2 * PI * to_int( f_AM / ( 2 * PI ) );
// calculate Kepler's equation for Earth
f_V = f_AM + 2 * f_ET * sin( f_AM ) + 1.25 * f_ET2 * sin( 2 * f_AM );
f_R = ( 1 - f_ET2 ) / ( 1 + f_ET * cos( f_V ) );
if( f_V < 0 ){ f_V = 2 * PI + f_V; }
f_L = f_P + f_V;
f_L -= 2 * PI * to_int( f_L / ( 2 * PI ) );
// calculate right ascension & declination
f_Z = ( f_J2 - f_S ) / 365.2422;
f_OB = deg2rad( 23.452294 - ( 0.46845 * f_Z + 0.00000059 * f_Z * f_Z ) / 3600.0 );
f_DC = asin( sin( f_OB ) * sin( f_L ) );
f_AR = acos( cos( f_L ) / cos( f_DC ) );
if( f_L > PI ){ f_AR = 2 * PI - f_AR; }
f_Z = f_R * sin( f_OB ) * sin( f_L );
f_OB = rad2deg( f_OB );
f_L = rad2deg( f_L );
f_AR = f_AR * 12 / PI;
// convert declination from radians to degrees
f_DC = rad2deg( f_DC );
return ({ f_AR, f_DC, f_L, f_OB, f_R });
}
/* EOF */

Code: [Select]
/* ** HEADER for ../funs/astro/earth_moon.c **
2005 Apr 19 - T. Cook
functions returning moon information
moon phase elements from astro:moonphase at cpan.org
moon rise/set elements from www.xylem.f2s.com/kepler/
moon position from users.zoominternet.net/~matto
-----
2011-Jul-17 - T. Cook began altering for DeadSouls
2010-Jun-01 - T. Cook reformatted for notepad++
2010-Jan-02 - T. Cook modified for TUB
2005-Sep-28 - T. Cook tidied up
to do:
moon phase not completely accurate to other sources; fix
*/
#pragma strong_types

//------------------------------------------------------------------//
#include "/realms/ardneh/area/customdefs.h"
#include "../../inc/astro.h"
#include "../../inc/math.h"

inherit MY_FUNS "/math.c";
inherit MY_FUNS "/astro/chronfuns.c";
inherit MY_FUNS "/lang/longnumbers.c";
inherit MY_FUNS "/geo/timezone.c";
/*
Dependencies:
ARDNEH_FUN"math.c"
ARDNEH_FUN"astro/chronfuns.c"
ARDNEH_FUN"lang/longnumbers.c"
ARDNEH_FUN"geo/timezone.c"
*/
float calcKepler(float ma, float ea);
float calcMeanMoonPhase(int JulianDay, float k);
float calcTrueMoonPhase(float k, float phase);
float *MoonPhaseHunt(int JulianDay);
float *calcMoonPhase(int JulianDay);
string MoonPhaseName(int JulianDay);
float *calcMoonPos(int JD,float dayfrac);
/* info for calcKepler()
func: solve the equation of Kepler
uses: deg2rad()
args: float mean anomaly
float eccentric anomaly
rtrn: float relation between polar coords and an elapsed time
*/
float calcKepler(float ma, float ea){
float EPSILON=0.000001, delta, m=deg2rad(ma), e=m;
do{
delta=e-ea*sin(e)-m;
e-=delta/(1-ea*cos(e));
}while(abs(delta)>EPSILON);
return e;
}
/* info for calcMeanMoonPhase()
func: calculate the time of mean new moon for given date
uses: deg2rad()
args: int Julian day
float k
rtrn: float Julian day of mean new moon
*/
float calcMeanMoonPhase(int JulianDay, float k){
float t=(JulianDay-2415020.0)/36525,  // time from epoch J1900.0
t2=t*t,                       // square for frequent use
t3=t2*t;                      // cube for frequent use
return 2415020.75933+Synmonth*k
+0.0001178*t2
-0.000000155*t3
+0.00033*sin(deg2rad(166.56+132.87*t-0.009173*t2));
}
/* info for calcTrueMoonPhase()
func: calculate the true, corrected phase time of moon
uses: deg2rad()
args: float K
float phase selector (0.0, 0.25, 0.5, 0.75)
rtrn: float true, corrected phase time in Julian days
*/
float calcTrueMoonPhase(float k, float phase){
int apcor;
float t, t2, t3, pt, m, mprime, f;
apcor=0;
k+=phase;           // add phase to new moon time
t=k/1236.85;        // time in Julian centuries from epoch J1900.0
t2=t*t;             // square for frequent use
t3=t2*t;            // cube for frequent use
// mean time of phase
pt=2415020.75933
+Synmonth*k
+0.0001178*t2
-0.000000155*t3
+0.00033*sin(deg2rad(166.56+132.87*t-0.009173*t2));
// Sun's mean anomaly
m=359.2242
+29.10535608*k
-0.0000333*t2
-0.00000347*t3;
// moon's mean anomaly
mprime=306.0253
+385.81691806*k
+0.0107306*t2
+0.00001236*t3;
// moon's argument of latitude
f=21.2964
+390.67050646*k
-0.0016528*t2
-0.00000239*t3;
// corrections for new and full moon
if((phase<0.01)||(abs(phase-0.5)<0.01)){
pt+=(0.1734-0.000393*t)*sin(deg2rad(m))
+0.0021*sin(deg2rad(2*m))
-0.4068*sin(deg2rad(mprime))
+0.0161*sin(deg2rad(2*mprime))
-0.0004*sin(deg2rad(3*mprime))
+0.0104*sin(deg2rad(2*f))
-0.0051*sin(deg2rad(m+mprime))
-0.0074*sin(deg2rad(m-mprime))
+0.0004*sin(deg2rad(2*f+m))
-0.0004*sin(deg2rad(2*f-m))
-0.0006*sin(deg2rad(2*f+mprime))
+0.0010*sin(deg2rad(2*f-mprime))
+0.0005*sin(deg2rad(m+2*mprime));
apcor=1;
}
else if((abs(phase-0.25)<0.01||(abs(phase-0.75)<0.01))){
pt+=(0.1721-0.0004*t)*sin(deg2rad(m))
+ 0.0021 * sin(deg2rad(2 * m))
- 0.6280 * sin(deg2rad(mprime))
+ 0.0089 * sin(deg2rad(2 * mprime))
- 0.0004 * sin(deg2rad(3 * mprime))
+ 0.0079 * sin(deg2rad(2 * f))
- 0.0119 * sin(deg2rad(m + mprime))
- 0.0047 * sin(deg2rad(m - mprime))
+ 0.0003 * sin(deg2rad(2 * f + m))
- 0.0004 * sin(deg2rad(2 * f - m))
- 0.0006 * sin(deg2rad(2 * f + mprime))
+ 0.0021 * sin(deg2rad(2 * f - mprime))
+ 0.0003 * sin(deg2rad(m + 2 * mprime))
+ 0.0004 * sin(deg2rad(m - 2 * mprime))
- 0.0003 * sin(deg2rad(2 * m + mprime));
// first quarter correction.
if(phase<0.5){
pt+=0.0028-0.0004*cos(deg2rad(m))+0.0003*cos(deg2rad(mprime));
}
// last quarter correction.
else{
pt+=-0.0028+0.0004*cos(deg2rad(m))-0.0003*cos(deg2rad(mprime));
}
apcor=1;
}
if (!apcor){
write("Error in calculation of true moon phase.\n");return 0;
}
return pt;
}
/* info for *MoonPhaseHunt()
func: find Julian date of phases of the moon which surround the
given date.  Five phases are found, starting and ending
with the new moons which bound the current lunation
args: int Julian day
rtrn: float array Phases containing values as Julian day
*/
float *MoonPhaseHunt(int JulianDay){
float sdate, adate, k1, k2, nt1, nt2;
int yy,mo,dd;
string a,b,MM;
sdate=JulianDay+0.0;
//*Phases=allocate(5),
adate     = sdate - 45;
sscanf(JD2UTC(JulianDay,0.0),"%s, %d--%s--%d %s UTC",a,yy,MM,dd,b);
mo=monthNumber(MM,0);
//int yy  = to_int(calcDateFromJD(sdate)[7..]),
//int mo = monthNumber(calcDateFromJD(sdate)[3..5]);
//int dd   = to_int(calcDateFromJD(sdate)[0..1]);
k1        = to_float(to_int((yy+((mo-1)*(1.0/12.0))-1900)*12.3685));
nt1       = calcMeanMoonPhase(to_int(adate),k1);
adate     = nt1;
while (1){
adate += Synmonth;
k2     = k1 + 1;
nt2    = calcMeanMoonPhase(to_int(adate),k2);
if ((nt1 <= sdate) && (nt2 > sdate)){break;}
nt1 = nt2;
k1  = k2;
}
return ({
calcTrueMoonPhase(k1, 0.0),
calcTrueMoonPhase(k1, 0.25),
calcTrueMoonPhase(k1, 0.5),
calcTrueMoonPhase(k1, 0.75),
calcTrueMoonPhase(k2, 0.0)
});
//return Phases;
}
/* info for *calcMoonPhase()
func: calculate the phase of the moon
args: int Julian day
rtrn: float array
0 terminator phase angle as percentage of full circle
(range is -0.25 .. 0.75, 0.5 is full, 0.0 is new
1 illuminated fraction of moon
2 moon's age in days
3 distance of moon from centre of Earth
4 angular diametre subtended by moon as seen from centre of Earth
*/
float *calcMoonPhase(int JulianDay){
float pdate, pphase, mage, dist, angdia, sudist, suangdia,
Day,N,M,Ec,Lambdasun,ml,MM,MN,Ev,Ae,A3,MmP,mEc,A4,lP,V,lPP,
NP,y,x,Lambdamoon,BetaM,MoonAge,MoonPhase,MoonDist,
MoonDFrac,MoonAng,MoonPar,F,SunDist,SunAng,mpfrac;
pdate = JulianDay+0.0;
pphase;    // illuminated fraction
mage;      // age of moon in days
dist;      // distance in kilometres
angdia;    // angular diameter in degrees
sudist;    // distance to Sun
suangdia;  // sun's angular diameter
//== Calculation of the Sun's position. ==
Day = pdate - Epoch;                      // date within epoch
N   = fixangle((360 / 365.2422) * Day);   // mean anomaly of the Sun
M   = fixangle(N + Elonge - Elongp);      // convert from perigee
  // co-ord to epoch 1980.0
Ec  = calcKepler(M,Eccent);               // solve eq. of Kepler
Ec  = sqrt((1+Eccent)/(1-Eccent))*tan(Ec/2);
Ec  = 2*rad2deg(atan(Ec));                // true anomaly
Lambdasun = fixangle(Ec + Elongp);        // Sun geocent. ecl. long.
F   = ((1+Eccent*cos(deg2rad(Ec)))/
(1-Eccent*Eccent));                   // orbit. dist. factor
SunDist = Sunsmax/F;                      // distance to Sun in km
SunAng  = F*Sunangsiz;                    // Sun angular size (deg)
//== Calculation of the Moon's position. ==
ml  = fixangle(13.1763966*Day+Mmlong);    // Moon mean lon
MM  = fixangle(ml-0.1114041*Day-Mmlongp); // Moon mean anomaly
MN  = fixangle(Mlnode-0.0529539*Day);     // Moon asc node mean lon
Ev  = 1.2739*sin(deg2rad(2*(ml-Lambdasun)-MM));  // evection
Ae  = 0.1858*sin(deg2rad(M));             // annual equation
A3  = 0.37*sin(deg2rad(M));               // correction term
MmP = MM+Ev-Ae-A3;                        // corredted anomaly
mEc = 6.2886*sin(deg2rad(MmP));           // correction eq of centre
A4  = 0.214*sin(deg2rad(2*MmP));          // another correction term
lP  = ml+Ev+mEc-Ae+A4;                    // corrected longitude
V   = 0.6583*sin(deg2rad(2*(lP-Lambdasun)));  // variation
lPP = lP+V;                               // true longitude
NP  = MN-0.16*sin(deg2rad(M));            // corrected node lon.
y   = sin(deg2rad(lPP-NP))*cos(deg2rad(Minc));  // y incl coord
x   = cos(deg2rad(lPP-NP));               // x inclination coord
Lambdamoon = rad2deg(atanb(y,x))+NP;      // ecliptic longitude
BetaM = rad2deg(asin(sin(deg2rad(lPP-NP))*
sin(deg2rad(Minc))));                   // ecliptic lat.
//== Calculation of the phase of the Moon. ==
MoonAge   = lPP-Lambdasun;                // age of moon in degrees
MoonPhase = (1-cos(deg2rad(MoonAge)))/2;  // phase of the moon
MoonDist  = (Msmax*(1-Mecc*Mecc))/
(1+Mecc*cos(deg2rad(MmP+mEc)));// moon dist fm c. of E.
MoonDFrac = MoonDist/Msmax;
MoonAng   = Mangsiz/MoonDFrac;            // moon's angular diam
MoonPar   = Mparallax/MoonDFrac;          // moon's parallax
pphase    = MoonPhase;
mage      = Synmonth*(fixangle(MoonAge)/360.0);
dist      = MoonDist;
angdia    = MoonAng;
sudist    = SunDist;
suangdia  = SunAng;
mpfrac    = fixangle(MoonAge)/360.0;
//float MoonDataArray;
//MoonDataArray=allocate(7);
//MoonDataArray=({mpfrac,pphase,mage,dist,angdia,sudist,suangdia});
//return MoonDataArray;
return ({mpfrac,pphase,mage,dist,angdia,sudist,suangdia});
}
/* info for MoonPhaseName()
func: get the name of the moon's phase
args: int Julian day
rtrn: string phase name
*/
string MoonPhaseName(int JulianDay){
float MP=calcMoonPhase(JulianDay)[0];
if     (MP>=-1.00&&MP<-0.73){return "first quarter";}
else if(MP>=-0.73&&MP<-0.53){return "waxing gibbous";}
else if(MP>=-0.53&&MP<-0.47){return "full";}
else if(MP>=-0.47&&MP<-0.27){return "waning gibbous";}
else if(MP>=-0.27&&MP<-0.23){return "third quarter";}
else if(MP>=-0.23&&MP<-0.03){return "waning crescent";}
else if(MP>=-0.03&&MP< 0.03){return "new";}
else if(MP>= 0.03&&MP< 0.23){return "waxing crescent";}
else if(MP>= 0.23&&MP< 0.27){return "first quarter";}
else if(MP>= 0.27&&MP< 0.47){return "waxing gibbous";}
else if(MP>= 0.47&&MP< 0.53){return "full";}
else if(MP>= 0.53&&MP< 0.73){return "waning gibbous";}
else if(MP>= 0.73&&MP< 1.00){return "third quarter";}
else                        {return "phantom";}
}
/* info for *calcMoonPos()
func: calculate moon position
args: int Julian day
float day fraction
rtrn: float array ({RA, decl., latitude, longitude, parallax})
*/
float *calcMoonPos(int JD,float dayfrac){
int DD,MM,YY,HR,MN;
string a,b,mo;
float T, T2, T3, L1, M, M1, D, F, OM;
float S, EX, L, B, W1, W2, BT, P, P1, LM, SEM, R, Z, OB, DC, AR, G1, S1;
sscanf(JD2UTC(JD,dayfrac),"%s, %d-%s-%([0-9][0-9]) %([0-9][0-9]):%([0-9][0-9]):%s UTC",a,YY,mo,DD,HR,MN,b);
DD = to_int( DD );
HR = to_int( HR );
MN = to_int( MN );
//<!--  LUNA   -->
//<!-- -->
T =(JD+dayfrac-2415020.0)/36525.0;
T2=T*T;
T3=T2*T;
L1=270.434164+481267.8831*T-0.001133*T2+0.0000019*T3;
M =358.475833+ 35999.0498*T-0.00015 *T2-0.0000033*T3;
M1=296.104608+477198.8491*T+0.009192*T2+0.0000144*T3;
D =350.737486+445267.1142*T-0.001436*T2+0.0000019*T3;
F = 11.250889+483202.0251*T-0.003211*T2-0.0000003*T3;
OM=deg2rad(259.183275-  1934.142 *T+0.002078*T2+0.0000022*T3);
//<!-- -->
L1=L1 + 0.000233*sin(deg2rad(51.2+20.2*T));
S=0.003964*sin(deg2rad(346.56+132.87*T-0.0091731*T2));
L1=L1+S+0.001964*sin(OM);
M =  M -0.001778*sin(deg2rad(51.2+20.2*T));
M1=  M1+0.000817*sin(deg2rad(51.2+20.2*T));
M1=M1+S+0.002541*sin(OM);
D =   D+0.002011*sin(deg2rad(51.2+20.2*T));
D = D+S+0.001964*sin(OM);
F = F+S-0.024691*sin(OM);
F =   F-0.004328*sin(OM+deg2rad(275.05-2.3*T));
EX =  1-0.002495*T-0.00000752*T2;
OM=rad2deg(OM);
L1=fixangle(L1);
M =fixangle(M);
M1=fixangle(M1);
D =fixangle(D);
F =fixangle(F);
OM=fixangle(OM);
M =deg2rad(M);
M1=deg2rad(M1);
D =deg2rad(D);
F =deg2rad(F);
//<!--longitude-->
L   = L1 + 6.288750*sin(M1)        +
       1.274018*sin(2*D-M1)    +
       0.658309*sin(2*D)       +
       0.213616*sin(2*M1)      -
    EX*0.185596*sin(M)         -
       0.114336*sin(2*F)       +
       0.058793*sin(2*D-2*M1)  +
    EX*0.057212*sin(2*D-M-M1)  +
       0.053320*sin(2*D+M1)    +
    EX*0.045874*sin(2*D-M)     +
    EX*0.041024*sin(M1-M)      -
       0.034718*sin(D)         -
    EX*0.030465*sin(M+M1)      +
       0.015326*sin(2*D-2*F)   -
       0.012528*sin(2*F+M1)    -
       0.010980*sin(2*F-M1)    +
       0.010674*sin(4*D-M1)    +
       0.010034*sin(3*M1)      +
       0.008548*sin(4*D-2*M1)  -
    EX*0.007910*sin(M-M1+2*D)  -
    EX*0.006783*sin(2*D+M)     +
       0.005162*sin(M1-D)      +
    EX*0.005000*sin(M+D)       +
    EX*0.004049*sin(M1-M+2*D)  +
       0.003996*sin(2*M1+2*D)  +
       0.003862*sin(4*D)       +
       0.003665*sin(2*D-3*M1)  +
    EX*0.002695*sin(2*M1-M)    +
       0.002602*sin(M1-2*F-2*D)+
    EX*0.002396*sin(2*D-M-2*M1)-
       0.002349*sin(M1+D)      +
EX*EX*0.002249*sin(2*D-2*M)   -
    EX*0.002125*sin(2*M1+M)    -
EX*EX*0.002079*sin(2*M)       +
EX*EX*0.002059*sin(2*D-M1-2*M)-
       0.001773*sin(M1+2*D-2*F)+
    EX*0.001220*sin(4*D-M-M1)  -
       0.001110*sin(2*M1+2*F)  +
       0.000892*sin(M1-3*D)    -
    EX*0.000811*sin(M+M1+2*D)  +
    EX*0.000761*sin(4*D-M-2*M1)+
EX*EX*0.000717*sin(M1-2*M)    +
EX*EX*0.000704*sin(M1-2*M-2*D)+
    EX*0.000693*sin(M-2*M1+2*D)+
    EX*0.000598*sin(2*D-M-2*F) +
       0.000550*sin(M1+4*D)    +
       0.000538*sin(4*M1)      +
    EX*0.000521*sin(4*D-M)     +
       0.000486*sin(2*M1-D)    -
       0.001595*sin(2*F+2*D);
//<!--latitude-->
B  =       5.128189*sin(F)         +
       0.280606*sin(M1+F)      +
       0.277693*sin(M1-F)      +
       0.173238*sin(2*D-F)     +
       0.055413*sin(2*D+F-M1)  +
       0.046272*sin(2*D-F-M1)  +
       0.032573*sin(2*D+F)     +
       0.017198*sin(2*M1+F)    +
       0.009266999*sin(2*D+M1-F)  +
       0.008823*sin(2*M1-F)    +
    EX*0.008247*sin(2*D-M-F)   +
       0.004323*sin(2*D-F-2*M1)+
       0.004200*sin(2*D+F+M1)  +
    EX*0.003372*sin(F-M-2*D)   +
    EX*0.002472*sin(2*D+F-M-M1)+
    EX*0.002222*sin(2*D+F-M)   +
       0.002072*sin(2*D-F-M-M1)+
    EX*0.001877*sin(F-M+M1)    +
       0.001828*sin(4*D-F-M1)  -
    EX*0.001803*sin(F+M)       -
       0.001750*sin(3*F)       +
    EX*0.001570*sin(M1-M-F)    -
       0.001487*sin(F+D)       -
    EX*0.001481*sin(F+M+M1)    +
    EX*0.001417*sin(F-M-M1)    +
    EX*0.001350*sin(F-M)       +
       0.001330*sin(F-D)       +
       0.001106*sin(F+3*M1)    +
       0.001020*sin(4*D-F)     +
       0.000833*sin(F+4*D-M1)  +
       0.000781*sin(M1-3*F)    +
       0.000670*sin(F+4*D-2*M1)+
       0.000606*sin(2*D-3*F)   +
       0.000597*sin(2*D+2*M1-F)+
    EX*0.000492*sin(2*D+M1-M-F)+
       0.000450*sin(2*M1-F-2*D)+
       0.000439*sin(3*M1-F)    +
       0.000423*sin(F+2*D+2*M1)+
       0.000422*sin(2*D-F-3*M1)-
    EX*0.000367*sin(M+F+2*D-M1)-
    EX*0.000353*sin(M+F+2*D)   +
       0.000331*sin(F+4*D)     +
    EX*0.000317*sin(2*D+F-M+M1)+
EX*EX*0.000306*sin(2*D-2*M-F) -
       0.000283*sin(M1+3*F);
W1 =       0.0004664 * cos(deg2rad(OM));
W2 =       0.0000754 * cos(deg2rad(OM + 275.05 - 2.3 * T));
BT = B * (1 - W1 - W2);
//<!--parallax-->
P  = 0.950724+0.051818*cos(M1)        +
          0.009531*cos(2*D-M1)    +
          0.007843*cos(2*D)       +
          0.002824*cos(2*M1)      +
          0.000857*cos(2*D+M1)    +
       EX*0.000533*cos(2*D-M)     +
       EX*0.000401*cos(2*D-M-M1)  +
          0.000173*cos(3*M1)      +
          0.000167*cos(4*D-M1)    -
       EX*0.000111*cos(M)         +
          0.000103*cos(4*D-2*M1)  -
          0.000084*cos(2*M1-2*D)  -
       EX*0.000083*cos(2*D+M)     +
          0.000079*cos(2*D+2*M1)  +
          0.000072*cos(4*D)       +
       EX*0.000064*cos(2*D-M+M1)  -
       EX*0.000063*cos(2*D+M-M1)  +
       EX*0.000041*cos(M+D)       +
       EX*0.000035*cos(2*M1-M)    -
          0.000033*cos(3*M1-2*D)  -
          0.000030*cos(M1+D)      -
          0.000029*cos(2*F-2*D)   -
       EX*0.000029*cos(2*M1+M)    +
    EX*EX*0.000026*cos(2*D-2*M)   -
          0.000023*cos(2*F-2*D+M1)+
       EX*0.000019*cos(4*D-M-M1);
P1 = P * 60;
B  = deg2rad(BT);
LM = deg2rad(L);
//<!-- -->
SEM= 10800* asin(deg2rad(0.272488 * P))/PI;
R  = 10800/(PI*P1);
//<!-- -->
Z  = (JD - 2415020.5) / 365.2422;
OB = deg2rad(23.452294 -(0.46845 * Z+ 0.00000059 * Z*Z)/3600);
DC =asin(sin(B) * cos(OB) + cos(B) * sin(OB) * sin(LM));
AR =acos(cos(B) * cos(LM) / cos(DC));
if (LM > PI){AR = 2 * PI - AR;}
B  = rad2deg(B);
LM = rad2deg(LM);
AR = AR * 12 / PI;
//<!-- -->
DC = rad2deg(DC);
//<!-- -->
D = abs(DC);
if (DC>0){G1=floor(D);}else{G1=(-1)*floor(D);}
M1 = floor((D - to_int(D)) * 60);
S1 = ((D - to_int(D)) * 60 - M1) * 60;
if (DC<0){M1=-M1;S1=-S1;}
return ({AR,DC,LM,B,P1});
}
/* EOF */

Code: [Select]
/* ** HEADER for ../funs/astro/planets.c **
2005-Jul-12 - T. Cook
positions of planets
(leeched from http://users.zoominternet.net/~matto/Java)
-----
2011-Jul-17 - T. Cook began altering for DeadSouls
2010-Jun-01 - T. Cook reformatted for notepad++
2010-Jan-02 - T. Cook modified for TUB, streamlined
2005-Sep-28 - T. Cook tidied up
*/
#pragma strong_types

#include "/realms/ardneh/area/customdefs.h"
#include "../../inc/astro.h"
#include "../../inc/math.h"

inherit MY_FUNS "/math.c";
inherit MY_FUNS "/astro/chronfuns.c";
inherit MY_FUNS "/lang/longnumbers.c";
inherit MY_FUNS "/geo/timezone.c";

float *calcPlanetPos(int JD,float dayfrac);
/* info for calcPlanetPos()
func: calculate position of planets
args: int JD, float dayfrac
rtrn: float array ({planet})({RA,dec,dist})
*/
float *calcPlanetPos(int JD,float dayfrac){
int HR,mm,ss;float MN;string a,b;
float *MA, *EX, *IN, *NOD, *VNOD, *T, *MLP, *MVLP, *MAM, *N;
float J2, *vT, J1, F1, U, OB;
float AR,DC,V1,V2,R,R1,R2,*XYZ2;
int j,i;
float *pp, W, *P, *Q, *XYZ1, *XYZ;
mixed *PlanetArray;
// Semimajor axes (AUs)
MA=({
0.3870986,
0.7233316,
1.5236883,
5.202561,
9.554747,
19.21814,
30.10957,
39.78459
});
// Orbital eccentricity (unitless)
EX=({
0.20563,
0.0067826,
0.0933865,
0.0484658,
0.0556155,
0.0463232,
0.0090021,
0.25387
});
// Orbital inclination (degrees)
IN=({
7.0043579,
3.394435,
1.8498011,
1.3041819,
2.4893741,
0.7729895,
1.7716016,
17.137
});
// Logitude of ascending nodes (degrees)
NOD=({
48.094173,
76.4997524,
49.4032001,
100.2520175,
113.4888341,
73.8768642,
131.5660649,
109.941
});
// Variations of ascending nodes
VNOD=({
0.00325,
0.00247,
0.00211,
0.00276,
0.00239,
0.00137,
0.00301,
0.00305
});
// Epoch of reference (Julian date)
T=({
2444238.5,
2444238.5,
2444238.5,
2444238.5,
2444238.5,
2444238.5,
2444238.5,
2444240.5
});
// Argument of perihelion (degrees)
MLP=({
77.14442128,
131.2895792,
335.690816,
14.0095493,
92.665397,
172.7363288,
47.8672148,
222.972
});
// Variation of periheliom
MVLP=({
0.004265,
0.003847,
0.005038,
0.004412,
0.005371,
0.004063,
0.003899,
0.00083
});
// Mean anomaly at epoch of reference (degrees)
MAM=({
154.15309,
224.44394,
150.61701,
132.95682,
72.65684,
55.33453,
212.49069,
346.467
});
// Mean motion (orbital revolutions/day)
N=({
4.092339,
1.602130,
0.524033,
0.0830912,
0.0334596,
0.0117308,
0.0059819,
0.0039746
});
// Julian date
sscanf(JD2UTC(JD,dayfrac),"%s %([0-9][0-9]):%([0-9][0-9])%([0-9][0-9]) UTC",HR,mm,ss);
HR = to_int( HR );
mm = to_int( mm );
ss = to_int( ss );
MN=mm+ss/60.0;
J2=JD+to_float(HR)/24.0+MN/1440.0;
// assign values for Earth
vT=({
282.60402,  // PT
0.00471,    // VPT
0.016718,   // ET
0.9856,     // NT
155.9044    // MT
});
// Julian date
J1=2444400.5;
F1=2415020.5;
// obliquity of eccentric
U = (J2 - F1) / 365.2422;
OB=23.45229444-(0.46845*U+0.00000059*U*U)/3600;
OB=deg2rad(OB);
// Kepler's Equation for Earth
vT[4]=deg2rad(vT[4]+vT[3]*(J2-J1)-to_int(vT[4]/360)*360);
V2=vT[4]+2*vT[2]*sin(vT[4])+1.25*vT[2]*vT[2]*sin(2*vT[4])+
(13*sin(3*vT[4])-3*sin(vT[4]))*vT[2]*vT[2]*vT[2]/12;
R2=(1-vT[2]*vT[2])/(1+vT[2]*cos(V2));
// cartesian coordinates for Earth
vT[0]=deg2rad(vT[0]+vT[1]*(J2-J1)/100);
XYZ2=({
R2*cos(vT[0]+V2),
R2*cos(OB)*sin(vT[0]+V2),
R2*sin(OB)*sin(vT[0]+V2)
});
pp=allocate(10);
P=allocate(3);
Q=allocate(3);
XYZ1=allocate(3);
XYZ=allocate(3);
PlanetArray=allocate(8);
for(i=0;i<8;i++){PlanetArray[i]=allocate(3);}
for(j=0;j<8;j++){
// data for the planet
pp=({
MA[j],       // ap
EX[j],       // ep
IN[j],       // i
NOD[j],      // nd
VNOD[j],     // vnd
T[j],        // t0
MLP[j],      // lp
MVLP[j],     // vlp
MAM[j],      // mp
N[j]         // np
});
// P & Q
pp[3]+=pp[4]*((J2 - J1)/100);
pp[6]+=pp[7]*((J2 - J1)/100);
W=deg2rad(pp[6]-pp[3]);
pp[3]=deg2rad(pp[3]);
pp[2]=deg2rad(pp[2]);
P[0]=cos(pp[3])*cos(W)-cos(pp[2])*sin(pp[3])*sin(W);
P[1]=sin(pp[3])*cos(W)*cos(OB)+cos(pp[2])*cos(pp[3])*sin(W)*cos(OB)-sin(pp[2])*sin(W) *sin(OB);
P[2]=sin(pp[3])*cos(W)*sin(OB)+cos(pp[2])*cos(pp[3])*sin(W)*sin(OB)+sin(pp[2])*sin(W) *cos(OB);
Q[0]=-cos(pp[3])*sin(W)-cos(pp[2])*sin(pp[3])*cos(W);
Q[1]=-sin(pp[3])*sin(W)*cos(OB)+cos(pp[2])*cos(pp[3])*cos(W)*cos(OB)-sin(pp[2])*cos(W)*sin(OB);
Q[2]=cos(pp[3])*cos(W)*sin(OB)*cos(pp[2])-sin(pp[3])*sin(W)*sin(OB)+sin(pp[2])*cos(W)*cos(OB);
// Kepler's Equation for the planet
pp[8]=deg2rad(pp[8]+pp[9]*(J2-pp[5])-to_int(pp[8]/360)*360);
V1=pp[8]+2*pp[1]*sin(pp[8])+1.25*pp[1]*pp[1]*sin(2*pp[8])+
(13*sin(3*pp[8])-3*sin(pp[8]))*pp[1]*pp[1]*pp[1]/12;
R1=pp[0]*(1-pp[1]*pp[1])/(1+pp[1]*cos(V1));
// planetary cartesian coordinates
XYZ1=({
P[0] * R1 * cos(V1) + Q[0] * R1 * sin(V1),
P[1] * R1 * cos(V1) + Q[1] * R1 * sin(V1),
P[2] * R1 * cos(V1) + Q[2] * R1 * sin(V1)
});
// geocentric cartesian coordinates
XYZ=({
XYZ1[0] + XYZ2[0],
XYZ1[1] + XYZ2[1],
XYZ1[2] + XYZ2[2]
});
// calculate Right Ascension  & Declination
AR=atan(XYZ[1]/XYZ[0])*12/PI;
DC=rad2deg(atan(XYZ[2]/sqrt(XYZ[0]*XYZ[0]+XYZ[1]*XYZ[1])));
// distance from earth
R=sqrt(XYZ[0]*XYZ[0]+XYZ[1]*XYZ[1]+XYZ[2]*XYZ[2]);
if(XYZ[0]<0){AR+=12;}
if(AR<0){AR+=24;}
// assign values to result array
PlanetArray[j]=({AR,DC,R});
}
return PlanetArray;
}
/* EOF */

Offline z993126

  • BFF
  • ***
  • Posts: 128
    • View Profile
Re: Sky map and astronomical functions
« Reply #4 on: July 23, 2011, 01:29:12 AM »
Code: [Select]
/* ** HEADER for ../funs/lang/longnumbers.c **
2005-Apr-19 - T. Cook
----
2011-Jul-17 - T. Cook began altering for DeadSouls
2010-Jun-01 - T. Cook reformatted for notepad++
2010-Jan-04 - T. Cook added spell_number2()
2009-Feb-11 - T. Cook finished fixing decimal spelling
2009-Jan-15 - T. Cook began revision for tublib
2005-Sep-27 - T. Cook tidied up
2005-May-25 - T. Cook finished adding ordinal-spelling ability
*/
#pragma strong_types

inherit "/realms/ardneh/funs/extras.c";

string ordinal(int number);
string spell_number2(int nr);
varargs string spell_number(mixed nr,int useordinal);
string groupname(int groupnumber);
/* info for ordinal()
func: get a number's ordinal
args: int number
rtrn: string "st","nd","rd","th"
*/
string ordinal(int number){
string n=to_string(number);
switch(n[<1..<1]){
case "1":if(sizeof(n)>1&&n[<2..<1]=="11"){return "th";}else{return "st";}
case "2":if(sizeof(n)>1&&n[<2..<1]=="12"){return "th";}else{return "nd";}
case "3":if(sizeof(n)>1&&n[<2..<1]=="13"){return "th";}else{return "rd";}
default:return "th";
}
}
/* info for spell_number2()
func: minimalist version
args: int number
rtrn: string spelled number
*/
string spell_number2(int number){
string *onetonine, *tentonineteen, *tentoninety, *nrgroups;
string var, *groupdigits, spelled, x;
int groups,currentgroup,i,a;
onetonine=({
"","one","two","three","four","five","six","seven","eight",
"nine"
});
tentonineteen=({
"","","","","","","","","","","ten","eleven","twelve",
"thirteen","fourteen","fifteen","sixteen","seventeen",
"eighteen","nineteen"
});
tentoninety=({
"","ten","twenty","thirty","forty","fifty","sixty","seventy",
"eighty","ninety"
});
nrgroups=({
"","","thousand","million","billion","quadrillion",
"quintillion","sextillion","septillion","octillion",
"nonillion","decillion","undecillion","duodecillion",
"tredecillion","quatturodecillion","quindecillion",
"sexdecillion","septendecillion","octodecillion",
"novemdecillion","vigintillion"
});
var=to_string(number);
spelled="";
if(var=="0"){return "zero";}else{
if(var[0..0]=="-"){spelled+="negative ";var=to_string(abs(number));}
groups=1+(strlen(var)-1)/3;
groupdigits=allocate(groups);
for(i=strlen(var);i>0;i-=3){
currentgroup=1+(strlen(var)/3)-(i/3);
if(currentgroup==groups){x=var[0..i-1];}else{x=var[i-3..i-1];}
groupdigits[currentgroup-1]=x;
}
}
for(i=groups-1;i>-1;i--){
if(groupdigits[i]!="000"){
if(strlen(groupdigits[i])>2){
if(i+1<groups){spelled+=" ";}
a=to_int(groupdigits[i][<3..<3]);
spelled+=onetonine[a]+" hundred";
}
if(strlen(groupdigits[i])>0){
if(strlen(groupdigits[i])>2){spelled+=" ";}
a=to_int(groupdigits[i][<2..<1]);
if(a<20&&a>9){
spelled+=tentonineteen[a];
}else{
a=to_int(groupdigits[i][<2..<2]);
spelled+=tentoninety[a];
if(strlen(groupdigits[i])>0){
a=to_int(groupdigits[i][<1..<1]);
if(a>1){spelled+="-";}
spelled+=onetonine[a];
}
}
}
if(i>0){spelled+=" "+nrgroups[i+1];}
}
}
return spelled;
}
/* info for spell_number()
func: spell a number
args: mixed number
int useordinal
rtrn: string spelled number, rounding if useordinal and float
*/
varargs string spell_number(mixed nr,int useordinal){
int a,n,i,x,hasdecimalold,decimalplaceold,z,oldnumber2;
string numberarrayold, b,c, *longones, *longonesord, newnumber, spellednumber;
int firstnonzero,lastnonzero,firstnumberfound,lastnumberfound;
int hasdecimal,decimalplace;
string numberarray,nr2;
float oldnumber1;
int predecimalgroups;
string thisgroupnumbers;
int currentgroup;
int firstnonzeronew,firstnumberfoundnew;
int postdecimalgroups;
string numberarraynew;
string pluralnumber;
n=0;
x=0;
hasdecimalold=0;
z=2147483646;
b="";
longones=({"zero","one","two","three","four","five","six","seven","eight","nine"});
longonesord=({"zeroth","first","second","third","fourth","fifth","sixth","seventh","eighth","ninth"});

// strip non-digits, make sure number can be spelled correctly
if(stringp(nr)){
if(nr[0..0]=="-"){b="-";}
for(i=0;i<sizeof(nr);++i){
//nr=regreplace(nr,"[^0-9\.]","",1);
if(sscanf(nr[i..i],"%d",c)==1){nr2+=nr[i..i];}
}
nr=b+nr2;
}
else if(nr>z||nr<-z||sizeof(regexp(({to_string(nr)}),"e+"))){
if(nr<0){b="negative ";}
switch(useordinal){
case 1:b+"big numberth";break;
default:b+="big number";break;
}
return b;
}
else if(sizeof(regexp(({to_string(nr)}),"e-"))){
if(nr<0){b="negative ";}b+="small number";return b;
}

numberarrayold=to_string(nr);  // convert number to string
if(!sizeof(numberarrayold)){return "not a number";}  // error
spellednumber="";  // initialize spelled number

// check for negative number, spell and remove from number
if(numberarrayold[0..0]=="-"){spellednumber+="negative ";n=1;}
newnumber=numberarrayold[n..];

// check for decimal point and its location
a=sizeof(newnumber);
for(i=0;i<a;i++){if(newnumber[i..i]=="."){hasdecimalold=1;decimalplaceold=i;}}
if(!hasdecimalold){decimalplaceold=sizeof(newnumber);a=decimalplaceold;}
if(useordinal){ // if using an ordinal, round any floats
oldnumber1=to_float(newnumber);
oldnumber2=to_int(newnumber[0..a]);
if(hasdecimalold){
if(oldnumber1-to_float(oldnumber2)>=0.5){
newnumber[a-1..a-1]=to_string(oldnumber2+1)[a-1..a-1];
}
newnumber=newnumber[0..a-1];
}
}

// strip padding zeroes
firstnonzero=0;
lastnonzero=0;
firstnumberfound=0;
lastnumberfound=0;
for(i=0;i<a;i++){if(newnumber[i..i]!="0"&&!firstnumberfound){firstnonzero=i;firstnumberfound=1;}}  // first nonzero
if(useordinal==0){  // last nonzero (skip if we're using ordinals)
for(i=sizeof(newnumber)-1;i>decimalplaceold;i--){if(newnumber[i..i]!="0"&&!lastnumberfound){lastnonzero=i;lastnumberfound=1;}}
if(lastnonzero==0){lastnonzero=decimalplaceold;}
}
else{lastnonzero=a;}
// put the stripped number into a new array
if(!useordinal){numberarray=newnumber[firstnonzero..lastnonzero];}
else{numberarray=newnumber[firstnonzero..decimalplaceold];}
// check for decimal point & location again; stripping padding zeroes changed where it is
a=sizeof(numberarray);
hasdecimal=0;
for(i=0;i<a;i++){if(numberarray[i..i]=="."){hasdecimal=1;decimalplace=i;}}
if(!hasdecimal){decimalplace=sizeof(numberarray);}
// if we're using an ordinal
if(useordinal){
// if using an ordinal, round any floats
oldnumber1=to_float(numberarray);
oldnumber2=to_int(numberarray[0..a]);
if(hasdecimal){
if(oldnumber1-to_float(oldnumber2)>=0.5){newnumber[a-1..a-1]=to_string(oldnumber2+1)[a-1..a-1];}
numberarray=numberarray[0..a-1];
}
}

// count the number of digit groups before the decimal point
predecimalgroups=(decimalplace+2)/3;
thisgroupnumbers="000";
currentgroup=(decimalplace-i+2)/3;
for(i=0;i<decimalplace;i++){
// put the groups digits into a string for handling
currentgroup=(decimalplace-i+2)/3;
switch((decimalplace-i)%3){
case 0:
thisgroupnumbers[0..0]=numberarray[i..i];
if(sizeof(numberarray)>1){thisgroupnumbers[1..1]=numberarray[i+1..i+1];}
if(sizeof(numberarray)>2){thisgroupnumbers[2..2]=numberarray[i+2..i+2];}
break;
case 2:
thisgroupnumbers[1..1]=numberarray[i..i];
if(sizeof(numberarray)>1){thisgroupnumbers[2..2]=numberarray[i+1..i+1];}
break;
case 1:
thisgroupnumbers[2..2]=numberarray[i..i];
break;
default:break;
}
switch((decimalplace-i)%3){
case 0:  // hundreds place
if(i%3==0){
thisgroupnumbers="000";
// put the groups digits into a string for handling
currentgroup=(decimalplace-i+2)/3;
switch((decimalplace-i)%3){
case 0:thisgroupnumbers[0..0]=numberarray[i..i];break;
case 2:thisgroupnumbers[1..1]=numberarray[i..i];break;
case 1:thisgroupnumbers[2..2]=numberarray[i..i];break;
default:break;
}
}
if(i>0&&numberarray[i..i]!="0"){spellednumber+=" ";}
if(numberarray[i..i]!="0"){spellednumber+=longones[to_int(numberarray[i..i])]+" hundred";}
break;
case 2:  // tens place
if(i>0&&numberarray[i..i]!="0"){spellednumber+=" ";}
switch(thisgroupnumbers[1..1]){
case "1":
switch(thisgroupnumbers[2..2]){
case "0":if(currentgroup==1&&useordinal!=0){spellednumber+="tenth";}else{spellednumber+="ten";}break;
case "1":spellednumber+="eleven";   break;
case "2":if(currentgroup==1&&useordinal!=0){spellednumber+="twelfth";}else{spellednumber+="twelve";}break;
case "3":spellednumber+="thirteen"; break;
case "4":spellednumber+="fourteen"; break;
case "5":spellednumber+="fifteen";  break;
case "6":spellednumber+="sixteen";  break;
case "7":spellednumber+="seventeen";break;
case "8":spellednumber+="eighteen"; break;
case "9":spellednumber+="nineteen"; break;
default: break;
}
break;
case "2":switch(thisgroupnumbers[2..2]){case "0":if(currentgroup==1&&useordinal!=0){spellednumber+="twentieth";}else{spellednumber+="twenty";}break;default: spellednumber+="twenty-";break;}break;
case "3":switch(thisgroupnumbers[2..2]){case "0":if(currentgroup==1&&useordinal!=0){spellednumber+="thirtieth";}else{spellednumber+="thirty";}break;default: spellednumber+="thirty-";break;}break;
case "4":switch(thisgroupnumbers[2..2]){case "0":if(currentgroup==1&&useordinal!=0){spellednumber+="fortieth";}else{spellednumber+="forty";}break;default: spellednumber+="forty-";break;}break;
case "5":switch(thisgroupnumbers[2..2]){case "0":if(currentgroup==1&&useordinal!=0){spellednumber+="fiftieth";}else{spellednumber+="fifty";}break;default: spellednumber+="fifty-";break;}break;
case "6":switch(thisgroupnumbers[2..2]){case "0":if(currentgroup==1&&useordinal!=0){spellednumber+="sixtieth";}else{spellednumber+="sixty";}break;default: spellednumber+="sixty-";break;}break;
case "7":switch(thisgroupnumbers[2..2]){case "0":if(currentgroup==1&&useordinal!=0){spellednumber+="seventieth";}else{spellednumber+="seventy";}break;default: spellednumber+="seventy-";break;}break;
case "8":switch(thisgroupnumbers[2..2]){case "0":if(currentgroup==1&&useordinal!=0){spellednumber+="eightieth";}else{spellednumber+="eighty";}break;default: spellednumber+="eighty-";break;}break;
case "9":switch(thisgroupnumbers[2..2]){case "0":if(currentgroup==1&&useordinal!=0){spellednumber+="ninetieth";}else{spellednumber+="ninety";}break;default: spellednumber+="ninety-";break;}break;
default:break;
}
break;
case 1:  // ones place
// if the number is 10 through 19, it's been checked
if(i>=1){
switch(thisgroupnumbers[1..1]){
case "1":break;
default:
switch(thisgroupnumbers[2..2]){
case "0":break;
default:
if(i>0&&thisgroupnumbers[2..2]!="0"&&thisgroupnumbers[1..1]=="0"){spellednumber+=" ";}
if(currentgroup==1&&useordinal!=0){
switch(thisgroupnumbers[2..2]){
case "4":case "6":case "7":spellednumber+=longones[to_int(numberarray[i..i])];break;
default:spellednumber+=longonesord[to_int(numberarray[i..i])];break;
}
}
else{spellednumber+=longones[to_int(numberarray[i..i])];}
break;
}
break;
}
}
else{
switch(thisgroupnumbers[2..2]){
case "0":
if(decimalplace==1&&numberarray[0..0]=="0"){
if(useordinal!=0){
if(currentgroup==1){spellednumber+=longonesord[to_int(numberarray[i..i])];}
else{spellednumber+=longones[to_int(numberarray[i..i])];}
}
else{spellednumber+=longones[to_int(numberarray[i..i])];}
}
break;
default:
if(i>0&&thisgroupnumbers[2..2]!="0"){spellednumber+=" ";}
if(currentgroup==1){
if(useordinal!=0){
switch(thisgroupnumbers[2..2]){
case "4":case "6":case "7":spellednumber+=longones[to_int(numberarray[i..i])];break;
default:spellednumber+=longonesord[to_int(numberarray[i..i])];break;
}
}
else{spellednumber+=longones[to_int(numberarray[i..i])];}
}
else{spellednumber+=longones[to_int(numberarray[i..i])];}
break;
}
}
// add the digit group name after the one if group wasn't 000
if(thisgroupnumbers!="000"){
if(currentgroup!=1){
b=groupname(currentgroup);
if(b!=""){spellednumber+=" "+b;}
else{spellednumber+=" "+to_string(decimalplace-i)+"-plex";}
}
}
// add the ordinal to numbers not accounted for yet
if(useordinal!=0){
if(currentgroup==1){
switch(sizeof(numberarray)){
case 1:
switch(thisgroupnumbers[2..2]){
case "0":spellednumber+="";break;
case "1":spellednumber+="";break;
case "2":spellednumber+="";break;
case "3":spellednumber+="";break;
case "5":spellednumber+="";break;
case "8":spellednumber+="";break;
case "9":spellednumber+="";break;
default: spellednumber+="th";break;
}
break;
default:
if((thisgroupnumbers[1..1]!="0"&&thisgroupnumbers[2..2]=="0")||
   (thisgroupnumbers[1..1]!="1"&&thisgroupnumbers[2..2]=="1")||
   (thisgroupnumbers[1..1]!="1"&&thisgroupnumbers[2..2]=="2")||
   (thisgroupnumbers[1..1]!="1"&&thisgroupnumbers[2..2]=="3")||
   (thisgroupnumbers[1..1]!="1"&&thisgroupnumbers[2..2]=="5")||
   (thisgroupnumbers[1..1]!="1"&&thisgroupnumbers[2..2]=="8")||
   (thisgroupnumbers[1..1]!="1"&&thisgroupnumbers[2..2]=="9")||
   (thisgroupnumbers[1..1]=="1"&&thisgroupnumbers[2..2]=="2")){
spellednumber+="";
}
else{spellednumber+="th";}
break;
}
}
}
break;
}
}
// can we combine these functions?
// number has a decimal place
if(hasdecimal&&!useordinal){
// check if there are digits before the decimal place
if(spellednumber!="zero"&&spellednumber!=""){spellednumber+=" and ";}
else{spellednumber="";}
// build a new array of just postdecimal digits
firstnumberfoundnew=0;
for(i=decimalplace+1,a=sizeof(numberarray);i<a;i++){
if(numberarray[i..i]!="0"&&!firstnumberfoundnew){firstnonzeronew=i;firstnumberfoundnew=1;}
}
numberarraynew=numberarray[decimalplace+1..];
// count postdecimal digit groups
postdecimalgroups=(sizeof(numberarraynew)+2)/3;
for(i=0,a=sizeof(numberarraynew);i<a;i++){
switch((a-i)%3){
case 0:  // hundreds place
if(i>0&&numberarraynew[i..i]!="0"){spellednumber+=" ";}
if(numberarraynew[i..i]!="0"){spellednumber+=longones[to_int(numberarraynew[i..i])]+" hundred";}
currentgroup=(a-i+2)/3;
break;
case 2:  // tens place
if(i>0&&numberarraynew[i..i]!="0"){spellednumber+=" ";}
switch(numberarraynew[i..i]){
case "1":
switch(numberarraynew[i+1..i+1]){
case "0":spellednumber+="ten";      break;
case "1":spellednumber+="eleven";   break;
case "2":spellednumber+="twelve";   break;
case "3":spellednumber+="thirteen"; break;
case "4":spellednumber+="fourteen"; break;
case "5":spellednumber+="fifteen";  break;
case "6":spellednumber+="sixteen";  break;
case "7":spellednumber+="seventeen";break;
case "8":spellednumber+="eighteen"; break;
case "9":spellednumber+="nineteen"; break;
default:break;
}
break;
case "2":switch(numberarraynew[i+1..i+1]){case "0":spellednumber+="twenty"; break;default: spellednumber+="twenty-";break;}break;
case "3":switch(numberarraynew[i+1..i+1]){case "0":spellednumber+="thirty"; break;default: spellednumber+="thirty-";break;}break;
case "4":switch(numberarraynew[i+1..i+1]){case "0":spellednumber+="forty"; break;default: spellednumber+="forty-";break;}break;
case "5":switch(numberarraynew[i+1..i+1]){case "0":spellednumber+="fifty"; break;default: spellednumber+="fifty-";break;}break;
case "6":switch(numberarraynew[i+1..i+1]){case "0":spellednumber+="sixty"; break;default: spellednumber+="sixty-";break;}break;
case "7":switch(numberarraynew[i+1..i+1]){case "0":spellednumber+="seventy"; break;default: spellednumber+="seventy-";break;}break;
case "8":switch(numberarraynew[i+1..i+1]){case "0":spellednumber+="eighty"; break;default: spellednumber+="eighty-";break;}break;
case "9":switch(numberarraynew[i+1..i+1]){case "0":spellednumber+="ninety"; break;default: spellednumber+="ninety-";break;}break;
default:break;
}
currentgroup=(a-i+2)/3;
break;
case 1:  // ones place
switch(numberarraynew[i-1..i-1]){
case "1":break;
default:
switch(numberarraynew[i..i]){
case "0":break;
default:
spellednumber+=longones[to_int(numberarraynew[i..i])];
break;
}
break;
}
// add the digit group name
currentgroup=(a-i+2)/3;
if(currentgroup<postdecimalgroups){thisgroupnumbers=to_string(numberarraynew[i-2..i]);}
if(thisgroupnumbers!="000"){
if(currentgroup!=1){
b=groupname(currentgroup);
if(b!=""){spellednumber+=" "+b;}
else{spellednumber+=" "+to_string(a-i)+"-plex";}
}
break;
}
default:break;
}
}
// add the decimal nths string
a=sizeof(numberarray);
switch((a-decimalplace-1)%3){
case 1: spellednumber+=" ten";    break;
case 2: spellednumber+=" hundred";break;
default:spellednumber+=" ";       break;
}
currentgroup=((a-decimalplace)+2)/3;
if(currentgroup>1&&(a-decimalplace-1)%3!=0){spellednumber+="-";}
if(currentgroup!=1){
b=groupname(currentgroup);
if(b!=""&&currentgroup!=1){spellednumber+=b;}
else{spellednumber+=" "+to_string(a-i)+"-plex";}
}
spellednumber+="th";
pluralnumber="";
if(replace_string(numberarraynew,"0","")!="1"){spellednumber+="s";}
}
return spellednumber;
}
/* info for groupname()
func: get digit group name
args: int groupnumber
rtrn: string spelled number group name
*/
string groupname(int groupnumber){
switch(groupnumber){
case   2:return "thousand";
case   3:return "million";
case   4:return "billion";
case   5:return "trillion";
case   6:return "quadrillion";
case   7:return "quintillion";
case   8:return "sextillion";
case   9:return "septillion";
case  10:return "octillion";
case  11:return "nonillion";
case  12:return "decillion";
case  13:return "undecillion";
case  14:return "duodecillion";
case  15:return "tredecillion";
case  16:return "quatturodecillion";
case  17:return "quindecillion";
case  18:return "sexdecillion";
case  19:return "septendecillion";
case  20:return "octodecillion";
case  21:return "novemdecillion";
case  22:return "vigintillion";
case 101:return "centillion";
default:return "";
}
}
/* EOF */


Offline daelaskai

  • BFF
  • ***
  • Posts: 174
    • View Profile
Re: Sky map and astronomical functions
« Reply #5 on: July 23, 2011, 11:29:16 AM »
Looks pretty cool.  Thanks for the contribution.

Daelaskai

Offline daelaskai

  • BFF
  • ***
  • Posts: 174
    • View Profile
Re: Sky map and astronomical functions
« Reply #6 on: July 23, 2011, 11:49:07 AM »
Looks like you're missing your directions.c reference file as noted in the sky.c object.

Offline z993126

  • BFF
  • ***
  • Posts: 128
    • View Profile
Re: Sky map and astronomical functions
« Reply #7 on: July 24, 2011, 03:00:34 PM »
(Important note:  one or more bits of code are deprecated and not actually needed in skymap.c ...this is one of them.  Had used at one point but since removed references, forgot to zap the #include when cleaning up.  Still, can be useful for other things, namely if you're wanting to reference the location of a star or such in the description of something.)

Code: [Select]
/* ** HEADER for ../funs/directions.c **
2005-Jun-28 - T. Cook
convert horizontal coordinates to world directions
-----
2011-Jul-17 - T. Cook began altering for DeadSouls
2010-May-31 - T. Cook reformatted for notepad++
2009-Mar-15 - T. Cook modified for tublib, tidied up
2005-Sep-27 - T. Cook tidied up, renamed functions
*/
#pragma strong_types

/* info for az_dir()
func: convert azimuth to compass direction
args: mixed azimuth in degrees
rtrn: string compass direction
*/
string az_dir(mixed az){
  if(az>=000.00&&az<=011.25){return "north";}
  if(az> 011.25&&az<=033.75){return "north by northeast";}
  if(az> 033.75&&az<=056.25){return "northeast";}
  if(az> 056.25&&az<=078.75){return "east by northeast";}
  if(az> 078.75&&az<=101.25){return "east";}
  if(az> 101.25&&az<=123.75){return "east by southeast";}
  if(az> 123.75&&az<=146.25){return "southeast";}
  if(az> 146.25&&az<=168.75){return "south by southeast";}
  if(az> 168.75&&az<=191.25){return "south";}
  if(az> 191.25&&az<=213.75){return "south by southwest";}
  if(az> 213.75&&az<=236.25){return "southwest";}
  if(az> 236.25&&az<=258.75){return "west by southwest";}
  if(az> 258.75&&az<=281.25){return "west";}
  if(az> 281.25&&az<=303.75){return "west by northwest";}
  if(az> 303.75&&az<=326.25){return "northwest";}
  if(az> 326.25&&az<=348.75){return "north by northwest";}
  if(az> 348.75&&az<=360.00){return "north";}
  return "thataway";
}
/* info for alt_dir()
func: convert altitude to direction
args: mixed altitude in degrees
rtrn: string direction
*/
string alt_dir(mixed alt){
  if(alt> 80&&alt<=90){return "overhead";}
  if(alt> 70&&alt<=80){return "nearly overhead";}
  if(alt> 60&&alt<=70){return "very high in the sky";}
  if(alt> 50&&alt<=60){return "high in the sky";}
  if(alt> 40&&alt<=50){return "in the middle of the sky";}
  if(alt> 30&&alt<=40){return "near the middle of the sky";}
  if(alt> 20&&alt<=30){return "low in the sky";}
  if(alt> 10&&alt<=20){return "very low in the sky";}
  if(alt>= 0&&alt<=10){return "on the horizon";}
  if(alt         <  0){return "below the horizon";}
  return "over there";
}
/* EOF */

Offline z993126

  • BFF
  • ***
  • Posts: 128
    • View Profile
Re: Sky map and astronomical functions
« Reply #8 on: July 24, 2011, 03:07:04 PM »
...and here's an extra file with constellation data.
Code: [Select]
/* ** HEADER for ../funs/astro/constellations.c **
2005 Jul 14 - T. Cook
constellation data; abbreviation, name, RA, declination
-----
2011-Jul-17 - T. Cook began altering for DeadSouls
2010-May-31 - T. Cook reformatted for notepad++
2010-Jan-05 - T. Cook made match with star data, made function
2009-Mar-15 - T. Cook simplified into single #define
2005-Sep-28 - T. Cook tidied up
*/
#ifndef constellations
#define constellations

mixed *CONSTELL(){
return({
({ 01.250, 39.000,"And","Andromeda" }),
({ 10.166,-38.433,"Ant","Antlia" }),
({ 16.166,-76.283,"Aps","Apus" }),
({ 22.716,-10.183,"Aqr","Aquarius" }),
({ 19.666,-02.500,"Aql","Aquila" }),
({ 17.383,-53.583,"Ara","Ara" }),
({ 03.000, 20.000,"Ari","Aries" }),
({ 06.000, 41.733,"Aur","Auriga" }),
({ 14.733, 30.716,"Boo","Bootes" }),
({ 05.000,-40.000,"Cae","Caelum" }),
({ 05.766, 70.266,"Cam","Camelopardalis" }),
({ 08.683, 20.150,"Cnc","Cancer" }),
({ 13.166, 41.150,"CVn","Canes Venatici" }),
({ 06.866,-21.983,"CMa","Canis Major" }),
({ 08.000, 05.000,"CMi","Canis Minor" }),
({ 21.016,-20.233,"Cap","Capricornus" }),
({ 08.766,-59.883,"Car","Carina" }),
({ 01.016, 62.200,"Cas","Cassiopeia" }),
({ 13.133,-45.966,"Cen","Centaurus" }),
({ 22.516, 71.583,"Cep","Cepheus" }),
({ 01.416,-11.350,"Cet","Cetus" }),
({ 10.583,-79.750,"Cha","Chamaeleon" }),
({ 15.083,-59.016,"Cir","Circinus" }),
({ 05.766,-35.283,"Col","Columba" }),
({ 12.766, 21.833,"Com","Coma Berenices" }),
({ 18.633,-41.483,"CrA","Corona Australis" }),
({ 16.000, 30.000,"CrB","Corona Borealis" }),
({ 12.466,-13.333,"Crv","Corvus" }),
({ 11.000,-15.000,"Crt","Crater" }),
({ 12.450,-59.966,"Crx","Crux" }),
({ 20.450, 42.033,"Cyg","Cygnus" }),
({ 20.700, 13.816,"Del","Delphinus" }),
({ 05.233,-63.700,"Dor","Dorado" }),
({ 17.750, 62.516,"Dra","Draco" }),
({ 21.000, 10.000,"Equ","Equuleus" }),
({ 03.250,-45.816,"Eri","Eridanus" }),
({ 02.783,-31.633,"For","Fornax" }),
({ 07.183, 22.683,"Gem","Gemini" }),
({ 22.000,-45.000,"Gru","Grus" }),
({ 17.333, 29.983,"Her","Hercules" }),
({ 03.116,-52.800,"Hor","Horologium" }),
({ 12.866,-35.800,"Hya","Hydra" }),
({ 02.000,-75.000,"Hyi","Hydrus" }),
({ 21.000,-55.000,"Ind","Indus" }),
({ 22.466, 44.816,"Lac","Lacerta" }),
({ 10.666, 16.450,"Leo","Leo" }),
({ 10.300, 35.166,"LMi","Leo Minor" }),
({ 05.583,-19.316,"Lep","Lepus" }),
({ 15.000,-15.000,"Lib","Libra" }),
({ 15.400,-42.966,"Lup","Lupus" }),
({ 08.033, 45.316,"Lyn","Lynx" }),
({ 18.833, 36.816,"Lyr","Lyra" }),
({ 05.000,-80.000,"Men","Mensa" }),
({ 21.000,-35.000,"Mic","Microscopium" }),
({ 07.15, -05.733,"Mon","Monoceros" }),
({ 12.333,-70.333,"Mus","Musca" }),
({ 16.050,-52.016,"Nor","Norma" }),
({ 22.000,-85.000,"Oct","Octans" }),
({ 17.183,-04.233,"Oph","Ophiuchus" }),
({ 05.583, 04.583,"Ori","Orion" }),
({ 19.683,-66.133,"Pav","Pavo" }),
({ 22.750, 19.533,"Peg","Pegasus" }),
({ 03.716, 41.766,"Per","Perseus" }),
({ 01.000,-50.000,"Phe","Phoenix" }),
({ 06.000,-55.000,"Pic","Pictor" }),
({ 00.516, 11.083,"Psc","Pisces" }),
({ 22.000,-30.000,"PsA","Piscis Austrinus" }),
({ 07.566,-39.383,"Pup","Puppis" }),
({ 08.966,-31.033,"Pyx","Pyxis" }),
({ 04.000,-60.000,"Ret","Reticulum" }),
({ 19.900, 18.783,"Sge","Sagitta" }),
({ 19.783,-20.766,"Sgr","Sagittarius" }),
({ 00.500,-32.350,"Scl","Sculptor" }),
({ 16.983,-37.166,"Sco","Scorpius" }),
({ 18.666,-10.300,"Sct","Scutum" }),
({ 17.183,-04.233,"Ser","Serpens" }),
({ 10.266,-02.416,"Sex","Sextans" }),
({ 04.766, 18.866,"Tau","Taurus" }),
({ 19.000,-50.000,"Tel","Telescopium" }),
({ 02.116, 32.033,"Tri","Triangulum" }),
({ 16.083,-65.916,"TrA","Triangulum Australe" }),
({ 00.133,-64.966,"Tuc","Tucana" }),
({ 11.166, 55.383,"UMa","Ursa Major" }),
({ 15.000, 70.000,"UMi","Ursa Minor" }),
({ 09.683,-47.450,"Vel","Vela" }),
({ 13.216,-03.733,"Vir","Virgo" }),
({ 08.000,-70.000,"Vol","Volans" }),
({ 20.333, 25.066,"Vul","Vulpecula" })
});
}
#endif
/* EOF */

Offline dravster

  • Acquaintance
  • *
  • Posts: 11
  • Hippy coder
    • View Profile
    • Elenia
Re: Sky map and astronomical functions
« Reply #9 on: July 24, 2011, 10:41:49 PM »
Bugfix for longnumbers:  (line 143 for the version I have)

Code: [Select]
   else if(nr>z||nr<-z||sizeof(regexp(({to_string(nr)}),"e+"))){
      if(nr<0){b="negative ";}
      switch(useordinal){
         case 1:b+="big numberth";break;
         default:b+="big number";break;
      }
      return b;


That is, case1 should have...

Code: [Select]
b+="big numberth"
and not
Code: [Select]
b+"big numberth";
Also chronfuns.c, in criertime():

Code: [Select]
         switch( i_m ){
         case 45: s_longtime += "a quarter"; break;
         case 59: s_longtime + "a minute"; break;
         default: s_longtime += spell_number( 60 - i_m, 0 ) + " minutes"; break;
         }

Again a typo'ed +  where it should be += on   case 59.


Also for regional_divisions.c I've added my town and county:

   case "gb-hrt":
      switch(lower_case(city)) {
         case "watford": return ({ 51.6667, -0.4000, 69.0 });
         default: return ({ 51.0000, -0.5000, 50.0 });
         }




Also in chronfuns, pragmas on my mud gives ERR("*Trying to put array in string\n"). I fixed
this by changing the definitions of the global arrays months and days to be mixed * (rather than
 string *, as they are both arrays of arrays, not arrays of strings!)

But apart from these eventTurnOn() is now returning a nice skymap!

Nice work duder :D

Regards

Drav.

Offline dravster

  • Acquaintance
  • *
  • Posts: 11
  • Hippy coder
    • View Profile
    • Elenia
Re: Sky map and astronomical functions
« Reply #10 on: July 24, 2011, 11:06:23 PM »
For completeion, I've just gone out and had a check and it has the right phase, declination/ascension for the moon, Jupiter and my birth constellation Aquariaus... can't check much more than that as dawn is only 40 or so mins away and it's already quite bright out!  :)

\o/ wheeee now I can check astro charts without having to leave my mud :D

Offline z993126

  • BFF
  • ***
  • Posts: 128
    • View Profile
Re: Sky map and astronomical functions
« Reply #11 on: July 25, 2011, 06:25:10 AM »
Thanks for the heads up on the typos and the extra data!  I need to get around to adding more stars to stars.c as well as locations to regional_divisions.c.  As an extra bonus, here's my tester object that happens to include a data table for the planets' exact(ish) locations and an extended listing of almanac data; different time conversions, sunrise and set times, etc.  I should probably add functions for moonrise and set times, too...

Code: [Select]
/* ** HEADER for ../obj/tester.c **
2009-Jan-15 - T. Cook
-----
2011-Jul-17 - T. Cook began altering for DeadSouls
*/
#pragma strong_types

#include <lib.h>
inherit LIB_ITEM;
inherit LIB_ACTIVATE;

#include "../customdefs.h"

#include "../../inc/astro.h"
#include "../../inc/math.h"

inherit MY_FUNS "/math.c";
inherit MY_FUNS "/extras.c";
inherit MY_FUNS "/lang/longnumbers.c";
inherit MY_FUNS "/geo/timezone.c";
inherit MY_FUNS "/geo/regional_divisions.c";
inherit MY_FUNS "/directions.c";
inherit MY_FUNS "/astro/chronfuns.c";
inherit MY_FUNS "/boxdrawing.c";
inherit MY_FUNS "/astro/earth_sun.c";
inherit MY_FUNS "/astro/earth_moon.c";
inherit MY_FUNS "/astro/planets.c";

static void create(){
item::create();
SetKeyName( "tester" );
SetId( ({ "test device", "device", "tester" }) );
SetShort( "a test device" );
SetLong(
"This is a general object used to test code.  Various features can be activated on it:\n" +
"  'fit'      fit text into a confined space\n" +
"  'arrays'   see how arrays get defined\n" +
"  'almanac'  display almanac information\n" +
"  'math'     check some math\n" +
"  'planets'  calculate the positions of the planets"
);
SetMass( 10 );
SetBaseCost( 0 );
}

void init(){
item::init();
}

int eventTurnOn( mixed m_args, mixed m_foo ){
// declare variables
int i_count, i_d, i_fit_end, i_fit_end2, i_H, i_m, i_max, i_min, i_s, i_w, i_x, i_y;
float f_dst, f_HA, f_LST;
float *fa_AA, *fa_MP, *fa_MP2;
string s_D, s_M, s_st, s_temp, s_ut, s_w;
string *sa_Pi, *sa_temp;
mixed *ma_JD, *ma_JDnew, *ma_locus, *ma_PA, *ma_sea, *ma_sJD, *ma_SP, *ma_temp;
//-----
// set location for dependent functions
//ma_locus = ({ -5.0, "US-FL", "Seminole County", "Oviedo" });
//ma_locus = ({ 0.0, "xx-xx", "", "", 70.0, 0.0, 0.0 });
ma_locus = ({ -6.0, "US-IA", "Henry County", "Wayland" });
ma_locus += city_loci( ma_locus[1], ma_locus[2], ma_locus[3] );

i_w = this_player()->GetScreen()[0] - 1;
s_w = to_string( i_w );
s_ut = Utime();
s_st = systime( ma_locus[0], ma_locus[1][0..1] );
sscanf( s_st, "%s, %d-%s-%([0-9][0-9]) %([0-9][0-9]):%([0-9][0-9]):%([0-9][0-9])", s_D, i_y, s_M, i_d, i_H, i_m, i_s );
i_d = to_int( i_d );
i_H = to_int( i_H );
i_m = to_int( i_m );
i_s = to_int( i_s );

switch( m_args ){
//=========================================================================
// this section tests placing text into a space which requires truncation
// variables used:
//   int i_x, i_y, i_fit_end, i_fit_end2, i_count
//   string sa_temp
//   mixed array ma_temp
//-------------------------------------------------------------------------
case "fit":
sa_temp[0] = "test!";
ma_temp = allocate( 6, "......" );
sa_temp = allocate( 2 );
i_x = random( 4 ) + 1;
i_y = random( 4 ) + 1;
i_fit_end = sizeof( sa_temp[0] );
i_fit_end2 = i_fit_end / 2;
write(
"Putting '" + sa_temp[0] + "' (length " + i_fit_end + ") centred at row " + i_y + ", column " + i_x + ".\n"
);
for( i_count = 0; i_count < i_fit_end; ++i_count ){
while( i_x - i_fit_end2 + i_count < 1 ){ ++i_count; }
if( i_x - i_fit_end2 + i_count > sizeof( ma_temp[0] ) - 2 ){ break; }
ma_temp[i_y][i_x - i_fit_end2 + i_count] = sa_temp[0][i_count];
}
for( i_count = 0; i_count < 6; ++i_count ){
write( ma_temp[i_count] + "\n" );
}
write( "Fit complete." );
return 1;
//=========================================================================
// this section tests assigning values to arrays
// variables used:
//   mixed array ma_temp
//-------------------------------------------------------------------------
case "arrays":
write( "Testing array allocation & assignment." );
ma_temp = allocate( 4 );
ma_temp[0] = allocate( 4, "x" );
ma_temp[1] = allocate( 4, "x" );
ma_temp[2] = allocate( 4, "x" );
ma_temp[3] = allocate( 4, "x" );
write( ma_temp[0][0] + " " + ma_temp[0][1] );
ma_temp[0] = ({ "y", "z" });
write( ma_temp[0][0] + " " + ma_temp[0][1] );
write( sprintf( "%O", ma_temp[0] ) );
return 1;
//=========================================================================
// this section displays almanac data
// variables used:
//   int i_H, i_m, i_w
//   float f_dst
//   string s_ut, s_st
//   string array sa_temp
//   mixed array ma_locus, ma_JD, ma_sea, ma_sJD, ma_JDnew
//-------------------------------------------------------------------------
case "almanac":
//=====
write( sprintf( "%" + s_w + "'" + box( 0022, 0 ) + "'s", "" ) );
write( sprintf(
"Latitude:     %d degree%s %d minute%s %.1f second%s %s",
abs( dd2dms( ma_locus[4] )[0] ),
abs( dd2dms( ma_locus[4] )[0] ) == 1 ? "" : "s",
abs( dd2dms( ma_locus[4] )[1] ),
abs( dd2dms( ma_locus[4] )[1] ) == 1 ? "" : "s",
abs( dd2dms( ma_locus[4] )[2] ),
abs( dd2dms( ma_locus[4] )[2] ) == 1 ? "" : "s",
ma_locus[4] >= 0.0 ? "North" : "South"
) );
write( sprintf(
"Longitude:    %d degree%s %d minute%s %.1f second%s %s",
abs( dd2dms( ma_locus[5] )[0] ),
abs( dd2dms( ma_locus[5] )[0] ) == 1 ? "" : "s",
abs( dd2dms( ma_locus[5] )[1] ),
abs( dd2dms( ma_locus[5] )[1] ) == 1 ? "" : "s",
abs( dd2dms( ma_locus[5] )[2] ),
abs( dd2dms( ma_locus[5] )[2] ) == 1 ? "" : "s",
ma_locus[5] >= 0.0 ? "East" : "West"
) );
write( sprintf( "Altitude:     %+.1f metres above sea level", ma_locus[6] ) );
f_dst = DST( s_ut, ma_locus[0], ma_locus[1][0..1] );
sa_temp = timezone( ma_locus[0], ma_locus[1][0..1], f_dst );
write( sprintf( "Time Zone:    %s (%s; DST offset %+.1f hour%s)", sa_temp[0], sa_temp[1], f_dst, abs( f_dst ) == 1 ? "" : "s" ) );
//-----
write( sprintf( "%" + s_w + "'" + box( 0011, 0 ) + "'s", "" ) );
write( sprintf( "UNIX time:    %+010d seconds", time() ) );
ma_JD = UTC2JD( s_ut );
write( sprintf( "Julian date:  %d.%s ", ma_JD[0], to_string( ma_JD[1] )[2..] ) );
//ma_JD = Utime2JD();
//write( sprintf( "UNIX to JD:    %d.%s", ma_JD[0], to_string( ma_JD[1] )[2..] ) );
write( "UTC:          " + s_ut );
//write( sprintf( "JD to UTC:     %s", JD2UTC( ma_JD[0], ma_JD[1] ) ) );
write( "Local:        " + s_st );
//write( sprintf( "UTC to Local:  %s", UTC2Local( s_ut, ma_locus[0], ma_locus[1] ) ) );
write( sprintf( "Time AM/PM:   %d:%02d %s", iwrap( i_H, 1, 12 ), i_m, ampm( i_H ) ) );
write( indent_string( "Town crier:   " + criertime( sprintf( "%02d:%02d", i_H, i_m ) ), i_w, 0, 13 ) );
write( indent_string( "Formal date:  " + crierdate( s_st, 0 ), i_w, 0, 13 ) );
ma_sea = season( ma_JD[0], ma_JD[1], ma_locus[4] );
write( "Season:       " + ma_sea[0] );
ma_sJD = cut_float( ma_sea[1] );
write( "  Began:      " + JD2UTC( ma_sJD[0], ma_sJD[1] )[0..15] );
ma_sJD = cut_float( ma_sea[2] );
write( "  Ends:       " + JD2UTC( ma_sJD[0], ma_sJD[1] )[0..15] );
//-----
write( sprintf( "%" + s_w + "'" + box( 0011, 0 ) + "'s", "" ) );
write( "Sidereal:     " + JD2LST( ma_JD[0], ma_JD[1], ma_locus[5] ) + "\n" );
ma_JDnew = calcSunRSJD( ma_JD[0], ma_JD[1], ma_locus[4], ma_locus[5], "astro", 1 );
write( indent_string(
"Astro:        " +
JD2Local( ma_JDnew[0], ma_JDnew[1], ma_locus[0], ma_locus[1] )[17..] +
"  Solar illumination begins to affect visiblity of stars",
i_w, 0, 27
) );
ma_JDnew = calcSunRSJD( ma_JD[0], ma_JD[1], ma_locus[4], ma_locus[5], "nautical", 1 );
write( indent_string(
"Nautical:     " +
JD2Local( ma_JDnew[0], ma_JDnew[1], ma_locus[0], ma_locus[1] )[17..] +
"  Able to discern horizon aboard a ship",
i_w, 0, 27
) );
ma_JDnew = calcSunRSJD( ma_JD[0], ma_JD[1], ma_locus[4], ma_locus[5], "civil", 1 );
write( indent_string(
"Civil:        " +
JD2Local( ma_JDnew[0], ma_JDnew[1], ma_locus[0], ma_locus[1] )[17..] +
"  Able to read unaided by artificial light outdoors",
i_w, 0, 27
) );
ma_JDnew = calcSunRSJD( ma_JD[0], ma_JD[1], ma_locus[4], ma_locus[5], "rise", 1 );
write( indent_string(
"Sunrise:      " +
JD2Local( ma_JDnew[0], ma_JDnew[1], ma_locus[0], ma_locus[1] )[17..] +
"  Upper rim of sun crosses horizon",
i_w, 0, 27
) );
ma_JDnew = calcSolNoonJD( ma_JD[0], ma_locus[5] );
write(
"Solar noon:   " +
JD2Local( ma_JDnew[0], ma_JDnew[1], ma_locus[0], ma_locus[1] )[17..] +
"  Sun at apex"
);
ma_JDnew = calcSunRSJD( ma_JD[0], ma_JD[1], ma_locus[4], ma_locus[5], "set", -1 );
write( indent_string(
"Sunset:       " +
JD2Local( ma_JDnew[0], ma_JDnew[1], ma_locus[0], ma_locus[1] )[17..] +
"  Upper rim of sun crosses horizon",
i_w, 0, 27
) );
ma_JDnew = calcSunRSJD( ma_JD[0], ma_JD[1], ma_locus[4], ma_locus[5], "civil", -1 );
write( indent_string(
"Civil:        " +
JD2Local( ma_JDnew[0], ma_JDnew[1], ma_locus[0], ma_locus[1] )[17..] +
"  Unable to read unaided by artificial light outdoors",
i_w, 0, 27
) );
ma_JDnew = calcSunRSJD( ma_JD[0], ma_JD[1], ma_locus[4], ma_locus[5], "nautical", -1 );
write( indent_string(
"Nautical:     " +
JD2Local( ma_JDnew[0], ma_JDnew[1], ma_locus[0], ma_locus[1] )[17..] +
"  Unable to discern horizon aboard a ship",
i_w, 0, 27
) );
ma_JDnew = calcSunRSJD( ma_JD[0], ma_JD[1], ma_locus[4], ma_locus[5], "astro", -1 );
write( indent_string(
"Astro:        " +
JD2Local( ma_JDnew[0], ma_JDnew[1], ma_locus[0], ma_locus[1] )[17..] +
"  No solar illumination visible",
i_w, 0, 27
) );
//=====
write( sprintf( "%" + s_w + "'" + box( 0022, 0 ) + "'s", "" ) );
return 1;
//=========================================================================
// this section displays a table of the positions of the sun, moon, planets, and pluto
// variables used:
//   int i_count
//   float f_LST, f_HA
//   float array fa_AA, fa_MP, fa_MP2
//   string s_temp, s_w, s_ut
//   string array sa_Pi
//   mixed array ma_JD, ma_SP, ma_locus, ma_PA
//-------------------------------------------------------------------------
case "planets":
// display border & header
s_temp = ( sprintf( "%" + s_w + "'" + box( 0022, 0 ) + "'s\n", "" ) );
s_temp += "         %^BOLD%^Rt. As.   Declin.%^RESET%^ ";
s_temp += box( 1100, 0 ) + " %^BOLD%^Hr.Ang.%^RESET%^ " + box( 1100, 0 );
s_temp += "  %^BOLD%^Alt.    Azimuth%^RESET%^ " + box( 1100, 0 );
s_temp += " %^BOLD%^Dist. (AU)%^RESET%^";
write( s_temp );
// assign the time for the calculations
ma_JD = UTC2JD( s_ut );
// calcate & display sun position
ma_SP = calcSunPos( ma_JD[0], ma_JD[1]);//, ma_locus[0], ma_locus[1][0..1] );
f_LST = HHmmss2mm( sprintf( "%s", JD2LST( ma_JD[0], ma_JD[1], ma_locus[5] )[0..7] ) );
f_HA = calcHourAngle( f_LST / 60.0, ma_SP[0] );
fa_AA = Equatorial2Horizontal( ma_locus[4], ma_SP[1], f_HA );
write( sprintf(
"  %s  %' '7.4f  %+' '8.4f%s" + box( 1100, 0 ) + " %+' '7.3f " + box( 1100, 0 ) +
" %+' '7.3f%s %' '7.3f%s" + box( 1100, 0 ) + " %' '8.5f",
"Sol: ",
ma_SP[0],
ma_SP[1],
" ",
f_HA,
fa_AA[0],
" ",
fa_AA[1],
" ",
ma_SP[4]
) );
// calculate & display moon position
fa_MP = calcMoonPos( ma_JD[0], ma_JD[1] );
f_HA = calcHourAngle( f_LST / 60.0, fa_MP[0] );
fa_AA = Equatorial2Horizontal( ma_locus[4], fa_MP[1], f_HA );
fa_MP2 = calcMoonPhase( ma_JD[0] );
write( sprintf(
"  %s:  %' '7.4f  %+' '8.4f%s" + box( 1100, 0 ) + " %+' '7.3f " + box( 1100, 0 ) +
" %+' '7.3f%s %' '7.3f%s" + box( 1100, 0 ) + " %' '8.5f (%s, %' '7.3f%% illuminated)",
"Luna",
fa_MP[0],
fa_MP[1],
" ",
f_HA,
fa_AA[0],
" ",
fa_AA[1],
" ",
fa_MP2[3] / 149598000.0,
MoonPhaseName( ma_JD[0] ),
( fa_MP2[1] ) * 100
) );
// calculate & display planet positions
sa_Pi = ({
"Merc",
"Venu",
"Mars",
"Jupi",
"Satu",
"Uran",
"Nept",
"Plut",
});
ma_PA = calcPlanetPos( ma_JD[0], ma_JD[1] );
for( i_count = 0; i_count < 8; ++i_count ){
f_HA = calcHourAngle( f_LST / 60.0, ma_PA[i_count][0] );
fa_AA = Equatorial2Horizontal( ma_locus[4], ma_PA[i_count][1], f_HA );
write( sprintf(
"  %4s:  %' '7.4f  %+' '8.4f%s" + box( 1100, 0 ) + " %+' '7.3f " + box( 1100, 0 ) +
" %+' '7.3f%s %' '7.3f%s" + box( 1100, 0 ) + " %' '8.5f",
sa_Pi[i_count],
ma_PA[i_count][0],
ma_PA[i_count][1],
" ",
f_HA,
fa_AA[0],
" ",
fa_AA[1],
" ",
ma_PA[i_count][2]
) );
}
// display border
write( sprintf( "%" + s_w + "'" + box( 0022, 0 ) + "'s\n", "" ) );
return 1;
//=========================================================================
// this section does some number things
// variables used:
//   int i_min, i_max, i_count
//   string s_temp
//-------------------------------------------------------------------------
case "math":
write( "Checking integer-wrap function:\n" );
i_min = 1;
i_max = 12;
s_temp = "";
for( i_count = -48; i_count < 25; ++i_count ){
s_temp = s_temp + " " + sprintf( "%+3d=%+3d  ", i_count, iwrap( i_count, i_min, i_max ) );
if( i_count % i_max == -i_min || i_count % i_max == i_max - 1 ){ s_temp += "\n"; }
}
write( s_temp + "\n" );
write( "Spelling some numbers:\n" );
i_count = 1;      write( "     " + i_count + "=" + spell_number2( i_count ) );
i_count = 12;     write( "    " + i_count + "=" + spell_number2( i_count ) );
i_count = 123;    write( "   " + i_count + "=" + spell_number2( i_count ) );
i_count = 1234;   write( "  " + i_count + "=" + spell_number2( i_count ) );
i_count = 12345;  write( " " + i_count + "=" + spell_number2( i_count ) );
i_count = 123456; write( i_count + "=" + spell_number2( i_count ) );
write( "    1.1  = " +spell_number( 1.1, 0 ) );
write( "    2.56 = " + spell_number( 2.56, 0 ) );
write( "    1st  = " + spell_number( 1.1, 1 ) + "\n" );
write( "Converting Julian dates to UTC at year & date boundary:\n" );
write( JD2UTC( 2455197, 0.49998 ) );
write( JD2UTC( 2455197, 0.49999 ) );
write( JD2UTC( 2455197, 0.50000 ) );
write( JD2UTC( 2455197, 0.50001 ) );
return 1;
//=========================================================================
// default case for nothing happening
//-------------------------------------------------------------------------
default:
return 0;
}
return 1;
}
/* EOF */

Offline daelaskai

  • BFF
  • ***
  • Posts: 174
    • View Profile
Re: Sky map and astronomical functions
« Reply #12 on: July 26, 2011, 08:31:39 PM »
It would be neat to have a modification of this and have a sky map with our own unique constellations corresponding to our individual MUD.  Great work.

Offline z993126

  • BFF
  • ***
  • Posts: 128
    • View Profile
Re: Sky map and astronomical functions
« Reply #13 on: July 27, 2011, 02:51:01 AM »
That's easy--just a matter of modifying stars.c and constellations.c.

Changing the planets is simply a matter of changing the values in planets.c for the orbital elements:
semimajor axis
orbital eccentricity
orbital inclination
longitude of ascending node
variation of ascending node
epoch of reference
argument of perihelion
variation of perihelion
mean anomaly at reference epoch
mean motion

The problem is I've not found data for the variations of ascending note and perihilion anywhere.  A possible solution is to replace the math in calcPlanetPos() with that at http://www.astro.uu.nl/~strous/AA/en/reken/hemelpositie.html or http://ssd.jpl.nasa.gov/txt/aprx_pos_planets.pdf (preferably the JPL approximations) which I hope to incorporate for comparison of efficiency and accuracy.

For the sun and moon...that's a lot more convoluted, and quite beyond me at this time.

Offline quixadhal

  • BFF
  • ***
  • Posts: 642
    • View Profile
    • WileyMUD
Re: Sky map and astronomical functions
« Reply #14 on: July 27, 2011, 02:49:48 PM »
How hard would it be to change the origin?

IE:  Rather than having the starmap centered on a view from Earth, what about a view from a star a few hundred light years away?  Obviously, your sun/moon/planet stuff wouldn't apply as you'd have different ones. :)