function rad(x) { return x*Math.PI/180 }
function grd(x) { return x*180/Math.PI }
function reduc(x) { return Math.PI*2*(x-Math.floor(x)) }

/////////////////////////////////////////////////

function Mond(JJJJ,MM,DD,hh,mm,ss,geografischeLaenge,geografischeBreite)
{

var deltaT=deltaTZeit(JJJJ); 
var t=JulDatum(JJJJ,MM,DD,hh,mm,ss,deltaT);
var T=t[0];
var T0=t[1];

var nut=nutation(T);
var NutationLaenge=nut[0];
var NutationEkliptik=nut[1];

var geogrBreite=geografischeBreite*Math.PI/180;
var geogrLaenge=geografischeLaenge*Math.PI/180;

var Sternzeit=sternzeit(T0,hh,mm,ss,geografischeLaenge); 

var m=mond_exakt(T,Sternzeit,geogrBreite,NutationLaenge,NutationEkliptik,deltaT);

return m
}

////////////////////////////////////////////////////////

function mond_exakt(T,SZ,geoBr,NutationLaenge,NutationEkliptik,deltaT) {
var m=new Array();

var m=mondkoordinaten(T,NutationLaenge,NutationEkliptik);
var ls=m[0];  // ekliptikale Länge
var bs=m[1];  // ekliptikale Breite
var ra=m[2];  // scheinbare Rektaszension
var de=m[3];  // scheinbare Deklination
var r=m[4];   // Abstand
var p=m[5];   // Paralaxe
var NutationLaenge=m[6];
var ekliptik=m[7];
var M=m[8];    // Konstante
var M1=m[9];   // Konstante
var D=m[10];   // Konstante
var F=m[11];   // Konstante
var E1=m[12];  // Konstante
var k=m[13];

var s=sonne(T,NutationLaenge); 
var Rerde=s[0];
var laengesonne=s[1];

var lib=physischeMonddaten(T,ra,de,r,Rerde,ls,bs,laengesonne,NutationLaenge,ekliptik,M,M1,D,F,E1);
var LibrationLaenge=lib[0]; // Libration in Länge
var LibrationBreite=lib[1]; // Libration in Breite
var SubLaenge=lib[2];       // Sublänge Sonne
var SubBreite=lib[3];       // Subbreite Sonne
var Colo= lib[4];           // Colongitude
var Posw= lib[5]; 

var lg=lichtgrenze(Colo);
   
var topo=mond_topo(geoBr,SZ,ra,de,r);
var ratopo=topo[0]; 
var detopo=topo[1];
var rtopo=topo[2];
var ptopo=Math.asin(6378.147/rtopo);

var tekl=ekliptikkoordinaten(ratopo,detopo,ekliptik);
var ltopo=tekl[0];
var btopo=tekl[1];

var libt=physischeMonddaten(T,ratopo,detopo,rtopo,Rerde,ltopo,btopo,laengesonne,NutationLaenge,ekliptik,M,M1,D,F,E1);
var LibrationLaengetopo=libt[0]; // topografische Libration in Länge
var LibrationBreitetopo=libt[1]; // topografische Libration in Breite
var Poswtopo=libt[5];
var Lichtgrenze=libt[6];

var hm=Math.asin(0.272493*Math.sin(p));
var hmt=Math.asin(0.272493*Math.sin(ptopo));

m[0]=ls; //eklipt Länge
m[1]=bs; //eklipt. Breite
m[2]=ra; //scheinb. Rektszension, geozentrisch
m[3]=de; //scheinb. Deklination,geozentrisch
m[4]=r; //Abstand geozentrisch
m[5]=p; // Parallaxe geozentrisch
m[6]=ratopo; //scheinbare Rektaszension topozentrisch
m[7]=detopo; //scheinbare Deklination topozentrisch
m[8]=rtopo; // Abstand topozentrisch
m[9]=ptopo; //Parallaxe topozentrisch
m[10]=LibrationLaenge; //Libration Länge geozentrische
m[11]=LibrationBreite; //Libration Breite geozentrisch
m[12]=SubLaenge; // Sublänge der Sonne
m[13]=SubBreite; // Subbreite der Sonne
m[14]=Colo; // Colongitude
m[15]=LibrationLaengetopo; //Libration Länge topozentrisch
m[16]=LibrationBreitetopo; //Libration Breite topozentrisch
m[17]=Posw;
m[18]=Poswtopo;
m[19]=lg; // Lichtgrenze
m[20]=k;
m[21]=SZ;
m[22]=T;
m[23]=ekliptik;
m[24]=hm;  // Halbmesser
m[25]=hmt; // Halbmesser topozentrisch
m[26]=laengesonne;
m[27]=deltaT;

return m 
}

/////////////////////////////////////////////////

function mondkoordinaten(T,NutationLaenge,NutationEkliptik)
{
var x=new Array();

var WR = Math.PI / 180;
var RW = 180 / Math.PI;
var VK = Math.PI * 2;

var bk=bahnkonst(T);
var D=bk[0];
var F=bk[1];
var M=bk[2];
var M1=bk[3];
var L1=bk[4];

var Phasenwinkel =180-180*D/Math.PI -6.289*Math.sin(M1)+2.100*Math.sin(M)-1.274*Math.sin(2*D-M1)-0.658*Math.sin(2*D)-0.214*Math.sin(2*M1)-0.110*Math.sin(D);

var k = (1+Math.cos(Phasenwinkel*Math.PI/180))/2;

var st=stoerungen(T,D,F,M,M1,L1);
var dl=st[0];
var dr=st[1];
var db=st[2];
var E1=st[3];

var scheinbLaenge = L1 + rad(dl/1000000) + NutationLaenge;
var scheinbBreite = rad(db/1000000);
var Abstand = 385000.560 + dr/1000;
var Parallaxe = Math.asin(6378.147/Abstand);

var ekliptik=mittlere_ekliptik(T) + NutationEkliptik*Math.PI/180; // wahre

var rade=aequatorkoordinaten(scheinbLaenge,scheinbBreite,ekliptik);
x[0]=scheinbLaenge;
x[1]=scheinbBreite;
x[2]=rade[0]; // Rektaszension
x[3]=rade[1]; // Deklination
x[4]=Abstand;
x[5]=Parallaxe; 
x[6]=NutationLaenge;
x[7]=ekliptik // wahre Ekliptik
x[8]=M;
x[9]=M1;
x[10]=D;
x[11]=F;
x[12]=E1;
x[13]=k;
return x
}

///////////////////////////////////////////////7

function bahnkonst(T) {
var bk=new Array();

var T2=T*T;
var T3=T2*T;
var T4=T3*T;

var L1=218.3164591+481267.88134236*T-0.0013268*T2+T3/538841-T4/65194000;
L1=reduc(L1/360);

var D=297.8502042+445267.1115168*T-0.0016300*T2+T3/545868-T4/113065000;
D=reduc(D/360);

var M=357.5291092+35999.0502909*T-0.0001536*T2+T3/24490000;
M=reduc(M/360);

var M1=134.9634114+477198.8676313*T+0.0089970*T2+T3/69699-T4/14712000;
M1=reduc(M1/360);

var F=93.2720993+483202.0175273*T-0.0034029*T2-T3/3526000+T4/863310000;
F=reduc(F/360);

bk[0]=D;
bk[1]=F;
bk[2]=M;
bk[3]=M1;
bk[4]=L1;

return bk
}

/////////////////////////////

function stoerungen(T,D,F,M,M1,L1) {
var st=new Array();
var VK = Math.PI * 2;

var A1=119.75+131.849*T;
A1=reduc(A1/360);

var A2=53.09+479264.290*T;
A2=reduc(A2/360);

var A3=313.45+481266.484*T;
A3=reduc(A3/360);

var E1=1-0.002516*T-0.0000075*T*T;
var E2=E1*E1;

with (Math) {

var dl=
+6288774*sin(M1)
+1274027*sin(2*D-M1)
 +658314*sin(2*D)
 +213618*sin(2*M1)
 -185116*sin(M)*E1
 -114332*sin(2*F)
  +58793*sin(2*D-2*M1)
  +57066*sin(2*D-M-M1)*E1
  +53322*sin(2*D+M1)
  +45758*sin(2*D-M)*E1
  -40923*sin(M-M1)*E1
  -34720*sin(D)
  -30383*sin(M+M1)*E1
  +15327*sin(2*D-2*F)
  -12528*sin(M1+2*F)
  +10980*sin(M1-2*F)
  +10675*sin(4*D-M1)
  +10034*sin(3*M1)
  +8548*sin(4*D-2*M1)
  -7888*sin(2*D+M-M1)*E1
  -6766*sin(2*D+M)*E1
  -5163*sin(D-M1)
  +4987*sin(D+M)*E1
  +4036*sin(2*D-M+M1)*E1
  +3994*sin(2*D+2*M1)
  +3861*sin(4*D)
  +3665*sin(2*D-3*M1)
  -2689*sin(M-2*M1)*E1
  -2602*sin(2*D-M1+2*F)
  +2390*sin(2*D-M-2*M1)*E1
  -2348*sin(D+M1)
  +2236*sin(2*D-2*M)*E2
  -2120*sin(M+2*M1)*E1
  -2069*sin(2*M)*E2
  +2048*sin(2*D-2*M-M1)*E2
  -1773*sin(2*D+M1-2*F)
  -1595*sin(2*D+2*F)
  +1215*sin(4*D-M-M1)*E1
  -1110*sin(2*M1+2*F)
   -892*sin(3*D-M1)
   -810*sin(2*D+M+M1)*E1
   +759*sin(4*D-M-2*M1)*E1
   -713*sin(2*M-M1)*E1
   -700*sin(2*D+2*M-M1)*E2
   +691*sin(2*D+M-2*M1)*E1
   +596*sin(2*D-M-2*F)*E1
   +549*sin(4*D+M1)
   +537*sin(4*M1)
   +520*sin(4*D-M)*E1
   -487*sin(D-2*M1)
   -399*sin(2*D+M-2*F)*E1
   -381*sin(2*M1-2*F)
   +351*sin(D+M+M1)*E1
   -340*sin(3*D-2*M1)
   +330*sin(4*D-3*M1)
   +327*sin(2*D-M+2*M1)*E1
   -323*sin(2*M+M1)*E1
   +299*sin(D+M-M1)*E1
   +294*sin(2*D+3*M1);

dl=dl+3958*sin(A1)+1962*sin(L1-F)+318*sin(A2);

var dr=
-20905355*cos(M1)
-3699111*cos(2*D-M1)
-2955968*cos(2*D)
-569925*cos(2*M1)
+48888*cos(M)*E1
-3149*cos(2*F)
+246158*cos(2*D-2*M1)
-152138*cos(2*D-M-M1)*E1
-170733*cos(2*D+M1)
-204586*cos(2*D-M)*E1
-129620*cos(M-M1)*E1
+108743*cos(D)
+104755*cos(M+M1)*E1
+10321*cos(2*D-2*F)
+79661*cos(M1-2*F)
-34782*cos(4*D-M1)
-23210*cos(3*M1)
-21636*cos(4*D-2*M1)
+24208*cos(2*D+M-M1)*E1
+30824*cos(2*D+M)*E1
-8379*cos(D-M1)
-16675*cos(D+M)*E1
-12831*cos(2*D-M+M1)*E1
-10445*cos(2*D+2*M1)
-11650*cos(4*D)
+14403*cos(2*D-3*M1)
-7003*cos(M-2*M1)*E1
+10056*cos(2*D-M-2*M1)*E1
+6322*cos(D+M1)
-9884*cos(2*D-2*M)*E2
+5751*cos(M+2*M1)*E1
-4950*cos(2*D-2*M-M1)*E2
+4130*cos(2*D+M1-2*F)
-3958*cos(4*D-M-M1)*E1
+3258*cos(3*D-M1)
+2616*cos(2*D+M+M1)*E1
-1897*cos(4*D-M-2*M1)*E1
-2117*cos(2*M-M1)*E2
+2354*cos(2*D+2*M-M1)*E2
-1423*cos(4*D+M1)
-1117*cos(4*M1)
-1571*cos(4*D-M)*E1
-1739*cos(D-2*M1)
-4421*cos(2*M1-2*F)
+1165*cos(2*M+M1)*E2
+8752*cos(2*D-M1-2*F);

var db=
+5128122*sin(F)
 +280602*sin(M1+F)
 +277693*sin(M1-F)
 +173237*sin(2*D-F)
  +55413*sin(2*D-M1+F)
  +46271*sin(2*D-M1-F)
  +32573*sin(2*D+F)
  +17198*sin(2*M1+F)
   +9266*sin(2*D+M1-F)
   +8822*sin(2*M1-F)
   +8216*sin(2*D-M-F)*E1
   +4324*sin(2*D-2*M1-F)
   +4200*sin(2*D+M1+F)
   -3359*sin(2*D+M-F)*E1
   +2463*sin(2*D-M-M1+F)*E1
   +2211*sin(2*D-M+F)*E1
   +2065*sin(2*D-M-M1-F)*E1
   -1870*sin(M-M1-F)*E1
   +1828*sin(4*D-M1-F)
   -1794*sin(M+F)*E1
   -1749*sin(3*F)
   -1565*sin(M-M1+F)*E1
   -1491*sin(D+F)
   -1475*sin(M+M1+F)*E1
   -1410*sin(M+M1-F)*E1
   -1344*sin(M-F)*E1
   -1335*sin(D-F)
   +1107*sin(3*M1+F)
   +1021*sin(4*D-F)
    +833*sin(4*D-M1+F)
    +777*sin(M1-3*F)
    +671*sin(4*D-2*M1+F)
    +607*sin(2*D-3*F)
    +596*sin(2*D+2*M1-F)
    +491*sin(2*D-M+M1-F)*E1
    -451*sin(2*D-2*M1+F)
    +439*sin(3*M1-F)
    +422*sin(2*D+2*M1+F)
    +421*sin(2*D-3*M1-F)
    -366*sin(2*D+M-M1+F)*E1
    -351*sin(2*D+M+F)*E1
    +331*sin(4*D+F)
    +315*sin(2*D-M+M1+F)*E1
    +302*sin(2*D-2*M-F)*E2
    -283*sin(M1+3*F)
    -229*sin(2*D+M+M1-F)*E1
    +223*sin(D+M-F)*E1
    +223*sin(D+M+F)*E1
    -220*sin(M-2*M1-F)*E1
    -220*sin(2*D+M-M1-F)*E1
    -185*sin(D+M1+F)
    +181*sin(2*D-M-2*M1-F)*E1
    -177*sin(M+2*M1+F)*E1
    +176*sin(4*D-2*M1-F)
    +166*sin(4*D-M-M1-F)*E1
    -164*sin(D+M1-F)
    +132*sin(4*D+M1-F)
    -119*sin(D-M1-F)
    +115*sin(4*D-M-F)*E1
    +107*sin(2*D-2*M+F)*E2;

db=db-2235*sin(L1)+382*sin(A3)+175*sin(A1-F)+175*sin(A1+F)+127*sin(L1-M1)-115*sin(L1+M1);
}
st[0]=dl; 
st[1]=dr;
st[2]=db;
st[3]=E1;
return st
}

///////////////////////////////////////////////

function mond_topo(geoBr,SZ,alphamond,deltamond,km) {
with (Math) {
var t=new Array();
  var erdr=6378.14-21.38*sin(geoBr)*sin(geoBr);
  var geozentrBr=geoBr-rad(0.1924*sin(2*geoBr));
  var xe=erdr*cos(geozentrBr)*cos(SZ);
  var ye=erdr*cos(geozentrBr)*sin(SZ);
  var ze=erdr*sin(geozentrBr)
  var xm=km*cos(deltamond)*cos(alphamond);
  var ym=km*cos(deltamond)*sin(alphamond);
  var zm=km*sin(deltamond);
  var xtopo=xm-xe;
  var ytopo=ym-ye;
  var ztopo=zm-ze;
  var kmtopo=sqrt(xtopo*xtopo + ytopo*ytopo + ztopo*ztopo);
  alphat=atan(ytopo/xtopo);
  if ( xtopo < 0) alphat=alphat+PI
  deltat=asin(ztopo/kmtopo);
  }
  t[0]=alphat;
  t[1]=deltat;
  t[2]=kmtopo;
  return t
  }
  

////////////////////////////////////////////

function Sonnenhoehe(SubSonneLaenge,SubSonneBreite,Objektlaenge,Objektbreite)
{
with (Math) {
var Colongitude = rad(90) - SubSonneLaenge;
var Sonnenhoehe = asin(Math.sin(SubSonneBreite)*sin(Objektbreite) + cos(SubSonneBreite)*cos(Objektbreite)*sin(Colongitude + Objektlaenge));
return Sonnenhoehe
}
}

/////////////////////////////////////////////////

function physischeMonddaten(T,RA,DE,Abstand,Rerde,scheinbL,scheinbB,laengesonne,NutationLaenge,ekliptik,M,M1,D,F,E1)
{
var lib=new Array();
var T2=T*T;
var T3=T2*T;
var T4=T3*T;

var I = rad(1.54242);

var Mondknoten = 125.0445550 - 1934.1361849*T +0.0020762*T2 + T3/467410 -T4/60616000;
Mondknoten=reduc(Mondknoten/360);

var x = optischeLibration(scheinbL,scheinbB,NutationLaenge,Mondknoten,I,F);
//  Wgeo=W*180/Math.PI;
var LibrationOptischLaenge = x[0];
var LibrationOptischBreite = x[1];
var Ae = x[2];

var kst=konstlib(T,M,M1,D,F,E1,Mondknoten);
var Rho = kst[0]; 
var Sigma = kst[1];
var Tau = kst[2];

var ph=libphys(Tau,Sigma,Rho,Ae,LibrationOptischBreite)
var LibrationPhysischLaenge=ph[0];
var LibrationPhysischBreite=ph[1];

var LibrationLaenge = LibrationOptischLaenge + LibrationPhysischLaenge;
var LibrationBreite = LibrationOptischBreite + LibrationPhysischBreite;

// Positionswinkel der Achse
var Poswinkel = poswachse(I,RA,Sigma,ekliptik,scheinbB,Mondknoten,NutationLaenge);


// selenografische Position der Sonne

var g=ekliptikkoordinaten(RA,DE,ekliptik);
var letopo=g[0];
var betopo=g[1];

var s=pos_sonne(Abstand,Rerde,laengesonne,letopo,betopo,NutationLaenge,Mondknoten,I,Tau,Sigma,Rho,F);
var SubLaenge=s[0];
var SubBreite=s[1];
var co=810-grd(SubLaenge);
co=reduc(co/360); 

var lg=lichtgrenze(co);

lib[0]=LibrationLaenge;
lib[1]=LibrationBreite;
lib[2]=SubLaenge;
lib[3]=SubBreite;
lib[4]=co;
lib[5]=Poswinkel;
lib[6]=lg // in Grad

return lib
}

//////////////////////////////////////////

function optischeLibration(scheinbL,scheinbB,NutationLaenge,Mondknoten,I,F)
{
var s = new Array();
with (Math) {
var W = scheinbL - NutationLaenge - Mondknoten;
var y = cos(W)*cos(scheinbB);
var x = sin(W)*cos(scheinbB)*cos(I) - sin(scheinbB)*sin(I);
var A = atan(x/y);
if (y<0) A = A + PI
var LibrationOL = A-F; 
var LibrationOB =asin(-sin(W)*cos(scheinbB)*sin(I) - sin(scheinbB)*cos(I));
}
s[0]=LibrationOL;
s[1]=LibrationOB;
s[2]=A;
s[3]=W;
return s
}

/////////////////////////////////////////////

function konstlib(T,M,M1,D,F,E1,Mondknoten) {

var k = new Array();
var K1 = rad(119.75 + 131.849 * T);
var K2 = rad(72.56 + 20.186 * T);

with (Math) {
var  Rho = 
- 0.02752 * cos(M1) 
- 0.02245 * sin(F) 
+ 0.00684 * cos(M1 - 2 * F)
- 0.00293 * cos(2 * F) 
- 0.00085 * cos(2 * F - 2 * D)
- 0.00054 * cos(M1 - 2 * D) 
- 0.0002 * sin(M1 + F)
- 0.0002 * cos(M1 + 2 * F) 
- 0.0002 * cos(M1 - F)
+ 0.00014 * cos(M1 + 2 * F - 2 * D);
Rho=rad(Rho);

var Sigma =
- 0.02816 * sin(M1) 
+ 0.02244 * cos(F) 
- 0.00682 * sin(M1 - 2 * F)
- 0.00279 * sin(2 * F) 
- 0.00083 * sin(2 * F - 2 * D)
+ 0.00069 * sin(M1 - 2 * D) 
+ 0.00040 * cos(M1 + F)
- 0.00025 * sin(2 * M1) 
- 0.00023 * sin(M1 + 2 * F)
+ 0.00020 * cos(M1 - F) 
+ 0.00019 * sin(M - F)
+ 0.00013 * sin(M1 + 2 * F - 2 * D) 
- 0.0001 * cos(M1 - 3 * F);
Sigma=rad(Sigma);

var  Tau = 
0.0252 * E1 * sin(M)                 
+ 0.00473 * sin(2 * M1 - 2 * F) 
- 0.00467 * sin(M1)
+ 0.00396 * sin(K1) 
+ 0.00276 * sin(2 * M1 - 2 * D)
+ 0.00196 * sin(Mondknoten)                 
- 0.00183 * cos(M1 - F)
+ 0.00115 * sin(M1 - 2 * D) 
- 0.00096 * sin(M1 - D)
+ 0.00046 * sin(2 * F - 2 * D) 
- 0.00039 * sin(M1 - F)
- 0.00032 * sin(M1 - M - 2 * D) 
+ 0.00027 * sin(2 * M1 - M - 2 * D)
+ 0.00023 * sin(K2) 
- 0.00014 * sin(2 * D)
+ 0.00014 * cos(2 * M1 - 2 * F) 
- 0.00012 * sin(M - 2 * F)
- 0.00012 * sin(2 * M1) 
+ 0.00011 * sin(2 * M1 - 2 * M - 2 * D);
Tau=rad(Tau);
}
k[0]=Rho;
k[1]=Sigma;
k[2]=Tau;
return k
}

////////////////////////////////////////

function poswachse(i,RA,Sigma,ekliptik,Breite,Mondknoten,NutationLaenge) {
with (Math) {
var V = Mondknoten + NutationLaenge + Sigma/sin(i); 
var X = sin(i+Sigma) * sin(V);
var Y = sin(i+Sigma) * cos(V) * cos(ekliptik) - cos(i+Sigma) * sin(ekliptik);
var omega = atan(X/Y);
if (Y<0) omega = omega + PI
var x = sqrt(X*X + Y*Y) * cos(RA-omega);
var y = cos(Breite);
pw = asin(x/y);
}
return pw
}

//////////////////////////////////////////

function libphys(Tau,Sigma,Rho,Ae,LibrationOptischBreite) {
var x=new Array();
with (Math) {
var LibrationPhysischLaenge = -Tau + (Rho*cos(Ae) + Sigma*sin(Ae))*tan(LibrationOptischBreite);
var LibrationPhysischBreite = Sigma*cos(Ae) - Rho*sin(Ae);
x[0]=LibrationPhysischLaenge;
x[1]=LibrationPhysischBreite;
}
return x
}

//////////////////////////////////////////////

function pos_sonne(Abstand,Rerde,laengesonne,scheinbLaenge,scheinbBreite,NutationLaenge,Mondknoten,I,Tau,Sigma,Rho,F) {
var s=new Array();

var IndexEntfernung = Abstand/(Rerde*149597870);
LaengeHeliozentrisch = laengesonne  + Math.PI + IndexEntfernung * Math.cos(scheinbBreite) * Math.sin(laengesonne-scheinbLaenge);
BreiteHeliozentrisch = IndexEntfernung * scheinbBreite;
var s = optischeLibration(LaengeHeliozentrisch,BreiteHeliozentrisch,NutationLaenge,Mondknoten,I,F);

// var W=s[3];
var A=s[2];
var LibrationOL=s[0];
var LibrationOB=s[1];
// var Whelio=W*180/Math.PI;  /// wo verwendet?
var As = A;
var LibL = LibrationOL;
var LibB = LibrationOB;

var ph=libphys(Tau,Sigma,Rho,As,LibB);
var LibrationPhysischLaenge=ph[0];
var LibrationPhysischBreite=ph[1];

s[0]=SubSonneLaenge = LibL + ph[0];
s[1]=SubSonneBreite = LibB + ph[1];
return s
}

function sonne(T,NutationLaenge) 
{
var T2=T*T;
var T3=T2*T;
var s=new Array();

var Ls = (280.46645 + 36000.76983*T + 0.0003032*T2)/360;
Ls = 2*Math.PI*(Ls-Math.floor(Ls));
var Ms = (357.52910 + 35999.05030*T - 0.0001559*T2)*Math.PI/180;
var es = 0.016708617 - 0.000042037*T - 0.0000001236*T3;
var Cs = (1.914600 - 0.004817*T - 0.000014*T2)*Math.sin(Ms)
  + (0.019993- 0.000101*T)*Math.sin(2*Ms) + 0.000290*Math.sin(3*Ms);
Cs = Cs*Math.PI/180
var vs = Ms + Cs;
var Rerde = (1.000001018*(1-es*es))/(1+es*Math.cos(vs));
var laengesonne = Ls + Cs;
laengesonne = laengesonne -0.00569*Math.PI/180 + NutationLaenge;
s[0]=Rerde;
s[1]=laengesonne;
return s
}

function aequatorkoordinaten(L, B, EKL)  {
// alpha=k[0], delta=k[1]
// alles in RAD
with (Math) {
var k=new Array();
var x = cos(B)*cos(L);
var y = cos(EKL)*cos(B)*sin(L)-sin(EKL)*sin(B);
k[0] = atan(y/x); if (x<0) k[0]  = k[0] +PI
k[1] = (sin(EKL)*cos(B)*sin(L))+(cos(EKL)*sin(B));
k[1] = asin(k[1]);
}
return k;
}


function ekliptikkoordinaten(ra,de,ekl) {
var n=new Array();
with (Math) {
var x=sin(ra)*cos(ekl)+tan(de)*sin(ekl);
var y=cos(ra);
var l=atan(x/y); if (y<0) l=l+PI
var b=sin(de)*cos(ekl)-cos(de)*sin(ekl)*sin(ra);
b=asin(b);
n[0]=l;
n[1]=b;
return n
}
}

function nutation(T)
{
var x=new Array();
var T2=T*T;
var T3=T2*T;
var T4=T3*T;
var VK = Math.PI * 2;
var L1=218.3164591+481267.88134236*T-0.0013268*T2+T3/538841-T4/65194000;
L1=L1/360;
L1=VK*(L1-Math.floor(L1));
var Knoten = 125.04452 - 1934.136261*T + 0.0020708*T2 +T3/450000;
Knoten = Knoten/360;
Knoten = VK* (Knoten-Math.floor(Knoten));
var LaengeSonne = (280.4665 + 36000.7698*T)/360;
LaengeSonne = LaengeSonne-Math.floor(LaengeSonne)*VK;
var NutationLaenge= -17.20*Math.sin(Knoten)- 1.32*Math.sin(2*LaengeSonne) - 0.23*Math.sin(2*L1) + 0.21*Math.sin(2*Knoten);
x[0] = Math.PI*NutationLaenge/(3600*180);
var NutationEkliptik = 9.20*Math.cos(Knoten) + 0.57*Math.cos(2*LaengeSonne) + 0.10*Math.cos(2*L1) - 0.09*Math.cos(2*Knoten);
x[1] = Math.PI*NutationEkliptik/(360*180);
return x
}


function sternzeit(T,hh,mm,ss,geoLaenge) {
var sz=100.46061837 +36000.770053608*T +0.0003879333*T*T - T*T*T/38710000;
var oz=1.0027379093*(1*hh+mm/60+ss/3600);
var sz1=1*sz+15*oz;
var sz2=sz1-geoLaenge; 
var szc=reduc(sz2/360); 
return szc
}


function JulDatum(Jahr,Monat,Tag,hh,mm,ss,deltaT) {
var z=new Array();
   JJJJ=Jahr;
   MM=Monat;
   DD=Tag;
   j=JJJJ;
   m=MM;
   if (m<=2) {j--; m=m+12;}
   B = Math.floor(j/400)-Math.floor(j/100);
   if (JJJJ<=1582 && MM<=10 && DD<=4) B=-2
   var JD0=Math.floor(365.25*j)+Math.floor(30.6001*(m+1))+B+1720996.5+DD;
   var T0=(JD0-2451545)/36525;
   s=1*ss+1*deltaT;
   d = JD0-2451545+hh/24+mm/1440+s/(1440*60);
   T = d/36525;
   z[0]=T;
   z[1]=T0;
   return z
   }
   
function mittlere_ekliptik(T)
{
var ekl = - 46.8150*T -0.00059*T*T + 0.001813*T*T*T;
var Ekliptik = rad(23.43929  + ekl/3600);
return Ekliptik
}

function lichtgrenze(colongitude) {
var c=colongitude*180/Math.PI;
var lg=360-c;
if (lg>180) lg=180-c
if (lg>90) lg=lg-180
lg=Math.round(lg*10)/10;
return lg
}

function deltaTZeit(Jahr) {
var dt=0;
if (Jahr>=1900) 
{ 
var T=(Jahr-2000)/100;
dt=((((-339.84*T-512.52)*T-160.22)*T)+92.23)*T+71.28;
dt=Math.round(1*dt); 
}
if (Jahr>1985) dt=delta_t(Jahr)
if (Jahr>=2008) dt=66

return dt
}


function delta_t(Jahr) {
var dt9=0;
var ddt=new Array(1986,54.9,1987,55.8,1988,56,1989,56,1990,56.9,1991,57,1992,58.3,1993,59,1994,60.00,1994,61,1996,61,1997,62.3,1998,63.0,1999,63.5,2000,63.8,2001,64.1,2002,64.3,2003,65.2,2004,65.8,2005,66.0,2006,66.0,2007,66);
for (var i=0;i<=ddt.length;i=i+2) 
{
if (ddt[i]==Jahr) dt9=ddt[i+1]
}
return dt9
}