//Global variables
var DST="off";
var MonthDays=new Array(12);
	MonthDays[1]=31;
	MonthDays[2]=28;
	MonthDays[3]=31;
	MonthDays[4]=30;
	MonthDays[5]=31;
	MonthDays[6]=30;
	MonthDays[7]=31;
	MonthDays[8]=31;
	MonthDays[9]=30;
	MonthDays[10]=31;
	MonthDays[11]=30;
	MonthDays[12]=31;

var NotCurrent="false"; //Flag for calculation time
var InDay;
var InMonth;
var InYear;		
var InHour;
var InMinute;

var pphase=0; //Illuminated fraction 
var mage=0; //Age of moon in days 
var dist=0; //Distance in kilometres 
var angdia=0; //Angular diameter in degrees 
var sudist=0; //Distance to Sun 
var suangdia=0; //Sun's angular diameter
var lunation=new Array(4); //5 found phases for current lunation
var yy; var mm; var dd; //for JYEAR function
var hh; var mmm; var ss; //for JHMS function
var tbuf; //for PHASETIME function

//Astronomical constants
var epoch=2444238.5; //1980 January 0.0
//Constants defining the Sun's apparent orbit
var elonge=278.833540; //Ecliptic longitude of the Sun at epoch 1980.0
var elongp=282.596403; //Ecliptic longitude of the Sun at perigee
var eccent=0.016718; //Eccentricity of Earth's orbit
var sunsmax=1.495985e8; //Semi-major axis of Earth's orbit, km
var sunangsiz=0.533128; //Sun's angular size, degrees, at semi-major axis distance
//Elements of the Moon's orbit, epoch 1980.0
var mmlong=64.975464; //Moon's mean longitude at the epoch
var mmlongp=349.383063; //Mean longitude of the perigee at the epoch
var mlnode=151.950429; // Mean longitude of the node at the epoch 
var minc=5.145396; //Inclination of the Moon's orbit 
var mecc=0.054900; //Eccentricity of the Moon's orbit 
var mangsiz=0.5181; //Moon's angular size at distance a from Earth 
var msmax=384401.0; //Semi-major axis of Moon's orbit in km 
var mparallax=0.9507; //Parallax at distance a from Earth 
var synmonth=29.53058868; //Synodic month (new Moon to new Moon) 
var lunatbase=2423436.0; //Base date for E. W. Brown's numbered series of lunations (1923 January 16) 
//Properties of the Earth  
var earthrad=6378.16; //Radius of Earth in kilometres 
//var PI=3.14159265358979323846;  //Assume not near black hole nor in Tennessee 
//Handy mathematical functions  

function sgn(x) //Extract sign
{
    var cn=x;
	var ns=0;
	if(cn<0)
		{ns=-1;}
	else
		{if(cn>0){ns=1;}}
	return ns;
}

function fixangle(a) //Fix angle
{
	var ca=a;
	var fa=0;
	fa=ca-360*Math.floor(ca/360);
	return fa;
}

function torad(d) //Deg->Rad
{
    var cd=d;
	var cr=0;
	cr=cd*(Math.PI/180);
	return cr;
}

function todeg(r) //Rad->Deg
{
	var cr=r;
	var cd=0;
	cd=cr*(180/Math.PI);
	return cd;
}

function dsin(x) //Sin from deg
{
	var cn=x;
	var nsin=0;
	nsin=Math.sin(torad(cn));
	return nsin;
}

function dcos(x) //Cos from deg
{
	var cn=x;
	var ncos=0;
	ncos=Math.cos(torad(cn));
	return ncos;
}

function julian() //Calculate Julian day
{	
	var CurYear=CurDate.getFullYear();
	var CurMonth=CurDate.getMonth()+1;
	var CurDay=CurDate.getDate();
	var CurHour=CurDate.getHours();
	var CurMinute=CurDate.getMinutes();
	var CurSecond=CurDate.getSeconds();
	var k1; var k2; var k3;
	var JulianDay;
	var HourCor;
	var JulianMidnight=12-(ActTZO/60);
	
	if(NotCurrent=="true"){
		JulianMidnight=14;
		CurDay=InDay;
		CurMonth=InMonth;
		CurYear=InYear;		
		CurHour=InHour;
		CurMinute=InMinute;
		CurSecond=0;
		if(document.CalcTime.IDST.status){DST="on";}
	}
	else{
		if(ActTZO==-180){DST="on"; JulianMidnight--;}
	}
				
	if(DST=="on"){CurHour=CurHour-1;}
	if(CurHour<0){
		CurDay=CurDay-1;
		if(CurDay<1){
			CurMonth=CurMonth-1;
			if(CurMonth<1){
				CurMonth=12;
				CurYear=CurYear-1;}
			CurDay=MonthDays[CurMonth];}
		CurHour=23;}
	
	CurHour=CurHour+(CurMinute/60)+(CurSecond/3600);

	CurYear=CurYear-Math.floor((12-CurMonth)/10);
	CurMonth=CurMonth+9; 
	if(CurMonth>=12){CurMonth=CurMonth-12;}

	k1=Math.floor(365.25*(CurYear+4712));
	k2=Math.floor(30.6*CurMonth+0.5);
	k3=Math.floor(Math.floor((CurYear/100)+49)*0.75)-38;

	JulianDay=k1+k2+CurDay+59; 
	if(JulianDay>2299160){JulianDay=JulianDay-k3;}
	if(CurHour<JulianMidnight){JulianDay=JulianDay-1;}
	if(CurHour-JulianMidnight<0)
		{HourCor=1-((1/24)*(JulianMidnight-CurHour));}
	else
		{HourCor=(1/24)*(CurHour-JulianMidnight);}

	JulianDay=JulianDay+HourCor;

	return JulianDay;
}

/*JYEAR - Convert	Julian	date  to  year,  month, day, which are
returned via integer pointers to integers (note that year is a long).*/
function jyear(tdi)
{
	var td=tdi+0.5;
	var z; var f; var a; var alpha;
	var b; var c; var d; var ee;

	z=Math.floor(td);
	f=td-z;
	if(z<2299161.0){a=z;} 
	else {
		alpha=Math.floor((z-1867216.25)/36524.25);
		a=z+1+alpha-Math.floor(alpha/4);
	}
	b=a+1524;
	c=Math.floor((b-122.1)/365.25);
	d=Math.floor(365.25*c);
	ee=Math.floor((b-d)/30.6001);
	dd=Math.floor(b-d-Math.floor(30.6001*ee)+f);
	if(ee<14){mm=ee-1;}
	else {mm=ee-13;}
	if(mm>2){yy=c-4716;}
	else {yy=c-4715;}
}

//JHMS - Convert Julian time to hour, minutes, and seconds.  
function jhms(j)
{
    var ij;
    j+=0.5; //Astronomical to civil 
    ij=(j-Math.floor(j))*86400.0+0.5;  //Round to nearest second
    hh=Math.floor(ij/3600)+Math.abs(ActTZO/60);
    mmm=Math.floor((ij/60)/60);
    ss=Math.floor(ij/60);
}

/*PHASETIME - Format  the  provided  julian  date  into  the
provided  buffer  in  UTC  format  for  screen display*/  
function phasetime(utime)
{
    jyear(utime);
    jhms(utime);
	if(hh.toString().length==1){hh="0"+hh;}
	if(mmm.toString().length==1){mmm="0"+mmm;}
	if(dd.toString().length==1){dd="0"+dd;}
	if(mm.toString().length==1){mm="0"+mm;}
	tbuf=dd+"."+mm+"."+yy+"   "+hh+":"+mmm;
}

/*MEANPHASE - Calculates  time  of  the mean new Moon for a given
base date.  This argument K to this function is the
precomputed synodic month index, given by:
K = (year - 1900) * 12.3685
where year is expressed as a year and fractional year.*/
function meanphase(sdate, k)
{
    var t; var t2; var t3; var nt11;

    //Time in Julian centuries from 1900 January 0.5 
    t=(sdate-2415020.0)/36525;
    t2=t*t; //Square for frequent use 
    t3=t2*t; //Cube for frequent use 
    nt11=2415020.75933+synmonth*k+0.0001178*t2-0.000000155*t3+0.00033*dsin(166.56+132.87*t-0.009173*t2);
    return nt11;
}

/*TRUEPHASE - Given a K value used to determine the mean phase of
the new moon, and a phase selector (0.0, 0.25, 0.5, 0.75), obtain 
the true, corrected phase time.*/
function truephase(k, phase)
{
    var t; var t2; var t3; var pt; var m; var mprime; var f;
    var apcor="FALSE";

	k+=phase; //Add phase to new moon time 
    t=k/1236.85; //Time in Julian centuries from 1900 January 0.5 
    t2=t*t; //Square for frequent use 
    t3=t2*t; //Cube for frequent use 
    pt=2415020.75933+synmonth*k+0.0001178*t2-0.000000155*t3+0.00033*dsin(166.56+132.87*t-0.009173*t2); //Mean time of phase 
    m=359.2242+29.10535608*k-0.0000333*t2-0.00000347*t3; //Sun's mean anomaly 
    mprime=306.0253+385.81691806*k+0.0107306*t2+0.00001236*t3; //Moon's mean anomaly 
    f=21.2964+390.67050646*k-0.0016528*t2-0.00000239*t3; //Moon's argument of latitude 
    if((phase<0.01) || (Math.abs(phase-0.5)<0.01)){
       //Corrections for New and Full Moon 
       pt+=(0.1734-0.000393*t)*dsin(m)+0.0021*dsin(2*m)-0.4068*dsin(mprime)+0.0161*dsin(2*mprime)-0.0004*dsin(3*mprime)+0.0104*dsin(2*f)-0.0051*dsin(m+mprime)-0.0074*dsin(m-mprime)+0.0004*dsin(2*f+m)-0.0004*dsin(2*f-m)-0.0006*dsin(2*f+mprime)+0.0010*dsin(2*f-mprime)+0.0005*dsin(m+2*mprime);
       apcor="TRUE";
    }
	else {
		if((Math.abs(phase-0.25)<0.01) || (Math.abs(phase-0.75)<0.01)){
		pt+=(0.1721-0.0004*t)*dsin(m)+0.0021*dsin(2*m)-0.6280*dsin(mprime)+0.0089*dsin(2*mprime)-0.0004*dsin(3*mprime)+0.0079*dsin(2*f)-0.0119*dsin(m+mprime)-0.0047*dsin(m-mprime)+0.0003*dsin(2*f+m)-0.0004*dsin(2*f-m)-0.0006*dsin(2*f+mprime)+0.0021*dsin(2*f-mprime)+0.0003*dsin(m+2*mprime)+0.0004*dsin(m-2*mprime)-0.0003*dsin(2*m+mprime);
		if(phase<0.5){pt+=0.0028-0.0004*dcos(m)+0.0003*dcos(mprime);} //First quarter correction 
		else {pt+=-0.0028+0.0004*dcos(m)-0.0003*dcos(mprime);} //Last quarter correction 
        apcor="TRUE";}
    }
    if(apcor="FALSE"){}
    return pt;
}

/*PHASEHUNT - Find time of phases of the moon which surround the
current date.  Five phases are found, starting and
ending with the new moons which bound the current
lunation.*/  
function phasehunt(sdate)
{
	var adate=sdate-45;
	var k1; var k2; var nt1; var nt2;
	var check="TRUE";

    jyear(adate);
    k1=Math.floor((yy+((mm-1)*(1/12))-1900)*12.3685);
	nt1=meanphase(adate, k1);
	adate=nt1;
    while(check="TRUE"){
        adate+=synmonth;
        k2=k1+1;
        nt2=meanphase(adate, k2);
        if (nt1<=sdate && nt2>sdate)
            break;
        nt1=nt2;
        k1=k2;}
    lunation[0] = truephase(k1, 0.0);
    lunation[1] = truephase(k1, 0.25);
    lunation[2] = truephase(k1, 0.5);
    lunation[3] = truephase(k1, 0.75);
    lunation[4] = truephase(k2, 0.0);
}

function kepler(im, iecc) //Solve the equation of Kepler.
{
    var m=im; var ecc=iecc;
	var ee; var delta; var EPSILON=1E-6;
    m=torad(m);
	ee=m; 
    do{
		delta=ee-ecc*Math.sin(ee)-m;
        ee-=delta/(1-ecc*Math.cos(ee));
    }while(Math.abs(delta)>EPSILON);
    return ee;
}

/*  PHASE - Calculate phase of moon as a fraction:
    The  argument  is  the  time  for  which  the  phase is requested,
    expressed as a Julian date and fraction.  Returns  the  terminator
    phase  angle  as a percentage of a full circle (i.e., 0 to 1), and
    stores into pointer arguments  the  illuminated  fraction  of  the
    Moon's  disc, the Moon's age in days and fraction, the distance of
    the Moon from the centre of the Earth, and  the  angular  diameter
    subtended  by the Moon as seen by an observer at the centre of the
    Earth.*/
function Phase(pdate) //pdate - Date for which to calculate phase 
{
    var Day; var N; var M; var Ec; var Lambdasun; var ml; var MM; var MN; var Ev;
	var Ae; var A3; var MmP; var mEc; var A4; var lP; var V; var lPP; var NP; 
	var y; var x; var Lambdamoon; var BetaM; var MoonAge; var MoonPhase; var MoonDist;
	var MoonDFrac; var MoonAng; var MoonPar; var F; var SunDist; var SunAng;

    //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-ordinates to epoch 1980.0 
    Ec=kepler(M, eccent); //Solve equation of Kepler 
    Ec=Math.sqrt((1+eccent)/(1-eccent))*Math.tan(Ec/2);
    Ec=2*todeg(Math.atan(Ec)); //True anomaly 
    Lambdasun=fixangle(Ec+elongp); //Sun's geocentric ecliptic longitude 
    //Orbital distance factor 
    F=((1+eccent*Math.cos(torad(Ec)))/(1-eccent*eccent));
    SunDist=sunsmax/F; //Distance to Sun in km 
    SunAng=F*sunangsiz; //Sun's angular size in degrees 
    //Calculation of the Moon's position 
    //Moon's mean longitude 
    ml=fixangle(13.1763966*Day+mmlong);
    //Moon's mean anomaly 
    MM=fixangle(ml-0.1114041*Day-mmlongp);
    //Moon's ascending node mean longitude 
    MN=fixangle(mlnode-0.0529539*Day);
    //Evection 
    Ev=1.2739*Math.sin(torad(2*(ml-Lambdasun)-MM));
    //Annual equation 
    Ae=0.1858*Math.sin(torad(M));
    //Correction term 
    A3=0.37*Math.sin(torad(M));
    //Corrected anomaly 
    MmP=MM+Ev-Ae-A3;
    //Correction for the equation of the centre 
    mEc=6.2886*Math.sin(torad(MmP));
    //Another correction term 
    A4=0.214*Math.sin(torad(2*MmP));
    //Corrected longitude 
    lP=ml+Ev+mEc-Ae+A4;
    //Variation 
    V=0.6583*Math.sin(torad(2*(lP - Lambdasun)));
    //True longitude 
    lPP=lP+V;
    //Corrected longitude of the node 
    NP=MN-0.16*Math.sin(torad(M));
    //Y inclination coordinate 
    y=Math.sin(torad(lPP-NP))*Math.cos(torad(minc));
    //X inclination coordinate 
    x=Math.cos(torad(lPP-NP));
    //Ecliptic longitude 
    Lambdamoon=todeg(Math.atan2(y,x));
    Lambdamoon+=NP;
    //Ecliptic latitude 
    BetaM=todeg(Math.asin(Math.sin(torad(lPP-NP))*Math.sin(torad(minc))));
    //Calculation of the phase of the Moon 
    //Age of the Moon in degrees 
    MoonAge=lPP-Lambdasun;
    //Phase of the Moon 
    MoonPhase=(1-Math.cos(torad(MoonAge)))/2;
    //Calculate distance of moon from the centre of the Earth 
    MoonDist=(msmax*(1-mecc*mecc))/(1+mecc*Math.cos(torad(MmP+mEc)));
    //Calculate Moon's angular diameter 
    MoonDFrac=MoonDist/msmax;
    MoonAng=mangsiz/MoonDFrac;
    //Calculate Moon's parallax 
    MoonPar=mparallax/MoonDFrac;

	pphase=MoonPhase;
    mage=synmonth*(fixangle(MoonAge)/360.0);
    dist=MoonDist;
    angdia=MoonAng;
    sudist=SunDist;
    suangdia=SunAng;
    return fixangle(MoonAge)/360.0;
}