From 79df146194e36e52e681aa859a06353e04139238 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Mon, 24 Jun 2019 14:08:38 -0400 Subject: [PATCH] More Python CalcMoon work in progress. Also made minor optimizations to spin() function in both JS and C. Fixed a mistake in JS and C that does not appear to have any algorithmic consequences, but it was definitely confusing once I saw it. --- generate/template/astronomy.c | 18 ++---- generate/template/astronomy.js | 2 +- generate/template/astronomy.py | 110 +++++++++++++++++++++++++++++++++ source/c/astronomy.c | 18 ++---- source/js/astronomy.js | 2 +- source/js/astronomy.min.js | 2 +- source/python/astronomy.py | 110 +++++++++++++++++++++++++++++++++ 7 files changed, 231 insertions(+), 31 deletions(-) diff --git a/generate/template/astronomy.c b/generate/template/astronomy.c index d5aa1c0b..8d381c55 100644 --- a/generate/template/astronomy.c +++ b/generate/template/astronomy.c @@ -1068,19 +1068,9 @@ static void spin(double angle, const double pos1[3], double vec2[3]) double angr = angle * DEG2RAD; double cosang = cos(angr); double sinang = sin(angr); - double xx = cosang; - double yx = sinang; - double zx = 0; - double xy = -sinang; - double yy = cosang; - double zy = 0; - double xz = 0; - double yz = 0; - double zz = 1; - - vec2[0] = xx*pos1[0] + yx*pos1[1] + zx*pos1[2]; - vec2[1] = xy*pos1[0] + yy*pos1[1] + zy*pos1[2]; - vec2[2] = xz*pos1[0] + yz*pos1[1] + zz*pos1[2]; + vec2[0] = +cosang*pos1[0] + sinang*pos1[1]; + vec2[1] = -sinang*pos1[0] + cosang*pos1[1]; + vec2[2] = pos1[2]; } static void ter2cel(astro_time_t time, const double vec1[3], double vec2[3]) @@ -1202,7 +1192,7 @@ static void Init(MoonContext *ctx) case 3: ARG=F; MAX=4; FAC=1.000002708+139.978*DGAM; break; case 4: ARG=D; MAX=6; FAC=1.0; break; } - CO(0,1) = 1.0; + CO(0,I) = 1.0; CO(1,I) = cos(ARG)*FAC; SI(0,I) = 0.0; SI(1,I) = sin(ARG)*FAC; diff --git a/generate/template/astronomy.js b/generate/template/astronomy.js index 20a3d65f..0babb994 100644 --- a/generate/template/astronomy.js +++ b/generate/template/astronomy.js @@ -752,7 +752,7 @@ function CalcMoon(time) { case 3: ARG=F; MAX=4; FAC=1.000002708+139.978*DGAM; break; case 4: ARG=D; MAX=6; FAC=1.0; break; } - SetCO(0, 1, 1); + SetCO(0, I, 1); SetCO(1, I, Math.cos(ARG) * FAC); SetSI(0, I, 0); SetSI(1, I, Math.sin(ARG) * FAC); diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index b3025540..7cd6f7e3 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -584,3 +584,113 @@ def _geo_pos(time, observer): pos2 = _nutation(time, -1, pos1) outpos = _precession(time.tt, pos2, 0.0) return outpos + +def _spin(angle, pos1): + angr = angle * _DEG2RAD + cosang = math.cos(angr) + sinang = math.sin(angr) + return [ + +cosang*pos1[0] + sinang*pos1[1], + -sinang*pos1[0] + cosang*pos1[1], + pos1[2] + ] + +def _ter2cel(time, vec1): + gast = _sidereal_time(time) + return _spin(-15.0 * gast, vec1) + +#---------------------------------------------------------------------------- +# CalcMoon + +class _Array1: + def __init__(xmin, xmax): + self.min = xmin + self.array = [0] * (xmax - xmin + 1) + + def __getitem__(self, key): + return self.array[key - self.min] + + def __setitem__(self, key, value): + self.array[key - self.min] = value + +class _Array2: + def __init__(xmin, xmax, ymin, ymax): + self.min = xmin + self.array = [_Array1(ymin, ymax) for i in range(xmax - xmin + 1)] + + def __getitem__(self, key): + return self.array[key - self.min] + + def __setitem__(self, key, value): + self.array[key - self.min] = value + +def _CalcMoon(time): + T = time.tt / 36525 + co = _Array2(-6, 6, 1, 4) + si = _Array2(-6, 6, 1, 4) + + def AddThe(c1, s1, c2, s2): + return (c1*c2 - s1*s2, s1*c2 + c1*s2) + + def Sine(phi): + return math.sin(_PI2 * phi) + + def Frac(x): + return x - math.floor(x) + + T2 = T*T + DLAM = 0 + DS = 0 + GAM1C = 0 + SINPI = 3422.7000 + S1 = Sine(0.19833+0.05611*T) + S2 = Sine(0.27869+0.04508*T) + S3 = Sine(0.16827-0.36903*T) + S4 = Sine(0.34734-5.37261*T) + S5 = Sine(0.10498-5.37899*T) + S6 = Sine(0.42681-0.41855*T) + S7 = Sine(0.14943-5.37511*T) + DL0 = 0.84*S1+0.31*S2+14.27*S3+ 7.26*S4+ 0.28*S5+0.24*S6 + DL = 2.94*S1+0.31*S2+14.27*S3+ 9.34*S4+ 1.12*S5+0.83*S6 + DLS =-6.40*S1 -1.89*S6 + DF = 0.21*S1+0.31*S2+14.27*S3-88.70*S4-15.30*S5+0.24*S6-1.86*S7 + DD = DL0-DLS + DGAM = ((-3332E-9 * Sine(0.59734-5.37261*T) + -539E-9 * Sine(0.35498-5.37899*T) + -64E-9 * Sine(0.39943-5.37511*T))) + + L0 = PI2*Frac(0.60643382+1336.85522467*T-0.00000313*T2) + DL0/ARC + L = PI2*Frac(0.37489701+1325.55240982*T+0.00002565*T2) + DL /ARC + LS = PI2*Frac(0.99312619+ 99.99735956*T-0.00000044*T2) + DLS/ARC + F = PI2*Frac(0.25909118+1342.22782980*T-0.00000892*T2) + DF /ARC + D = PI2*Frac(0.82736186+1236.85308708*T-0.00000397*T2) + DD /ARC + + I = 1 + while I <= 4: + if I == 1: + ARG=L; MAX=4; FAC=1.000002208 + elif I == 2: + ARG=LS; MAX=3; FAC=0.997504612-0.002495388*T + elif I == 3: + ARG=F; MAX=4; FAC=1.000002708+139.978*DGAM + else: + ARG=D; MAX=6; FAC=1.0 + + co[0][I] = 1 + co[1][I] = math.cos(ARG) * FAC + si[0][I] = 0 + si[1][I] = math.sin(ARG) * FAC + + J = 2 + while J <= MAX: + co[J][I], si[J][I] = AddThe(CO(J-1,I), SI(J-1,I), CO(1,I), SI(1,I)) + J += 1 + + J = 1 + while J <= MAX: + co[-J][I] = +co[J][I] + si[-J][I] = -si[J][I] + J += 1 + + I += 1 + diff --git a/source/c/astronomy.c b/source/c/astronomy.c index fe074f42..1f9773b8 100644 --- a/source/c/astronomy.c +++ b/source/c/astronomy.c @@ -1159,19 +1159,9 @@ static void spin(double angle, const double pos1[3], double vec2[3]) double angr = angle * DEG2RAD; double cosang = cos(angr); double sinang = sin(angr); - double xx = cosang; - double yx = sinang; - double zx = 0; - double xy = -sinang; - double yy = cosang; - double zy = 0; - double xz = 0; - double yz = 0; - double zz = 1; - - vec2[0] = xx*pos1[0] + yx*pos1[1] + zx*pos1[2]; - vec2[1] = xy*pos1[0] + yy*pos1[1] + zy*pos1[2]; - vec2[2] = xz*pos1[0] + yz*pos1[1] + zz*pos1[2]; + vec2[0] = +cosang*pos1[0] + sinang*pos1[1]; + vec2[1] = -sinang*pos1[0] + cosang*pos1[1]; + vec2[2] = pos1[2]; } static void ter2cel(astro_time_t time, const double vec1[3], double vec2[3]) @@ -1293,7 +1283,7 @@ static void Init(MoonContext *ctx) case 3: ARG=F; MAX=4; FAC=1.000002708+139.978*DGAM; break; case 4: ARG=D; MAX=6; FAC=1.0; break; } - CO(0,1) = 1.0; + CO(0,I) = 1.0; CO(1,I) = cos(ARG)*FAC; SI(0,I) = 0.0; SI(1,I) = sin(ARG)*FAC; diff --git a/source/js/astronomy.js b/source/js/astronomy.js index 474fbf4e..8f6d3105 100644 --- a/source/js/astronomy.js +++ b/source/js/astronomy.js @@ -1590,7 +1590,7 @@ function CalcMoon(time) { case 3: ARG=F; MAX=4; FAC=1.000002708+139.978*DGAM; break; case 4: ARG=D; MAX=6; FAC=1.0; break; } - SetCO(0, 1, 1); + SetCO(0, I, 1); SetCO(1, I, Math.cos(ARG) * FAC); SetSI(0, I, 0); SetSI(1, I, Math.sin(ARG) * FAC); diff --git a/source/js/astronomy.min.js b/source/js/astronomy.min.js index 0bd31137..ff6258ec 100644 --- a/source/js/astronomy.min.js +++ b/source/js/astronomy.min.js @@ -36,7 +36,7 @@ d+c>>1;if(br[g+1].mjd)d=g+1;else{b=r[g].dt+(b-r[g].mjd 4.848136811095359E-13*f)/4.84813681109536E-6;c=B(a);Q={tt:a.tt,dpsi:b,deps:f,ee:b*Math.cos(.017453292519943295*c)/15,mobl:c,tobl:c+f/3600}}return Q}function E(a){function b(a,b){var c=[],d;for(d=0;d<=b-a;++d)c.push(0);return{min:a,array:c}}function d(a,c,d,f){var g=[],e;for(e=0;e<=c-a;++e)g.push(b(d,f));return{min:a,array:g}}function c(a,b,c){a=a.array[b-a.min];return a.array[c-a.min]}function g(a,b,c){a=w.array[a-w.min];a.array[b-a.min]=c}function h(a,b,c){a=r.array[a-r.min];a.array[b-a.min]=c}function e(a, b,c,d,f){return f(a*c-b*d,b*c+a*d)}function l(a){return Math.sin(H*a)}function m(a,b,d,f){var g={x:1,y:0};a=[null,a,b,d,f];for(b=1;4>=b;++b)0!==a[b]&&e(g.x,g.y,c(w,a[b],b),c(r,a[b],b),function(a,b){return g.x=a,g.y=b});return g}function f(a,b,c,d,f,g,e,h){f=m(f,g,e,h);v+=a*f.y;q+=b*f.y;y+=c*f.x;A+=d*f.x}++I.calcmoon;a=a.tt/36525;var n,p,v,q,w=d(-6,6,1,4),r=d(-6,6,1,4);var u=a*a;var y=q=v=0;var A=3422.7;var z=l(.19833+.05611*a);var B=l(.27869+.04508*a);var K=l(.16827-.36903*a);var D=l(.34734-5.37261* a),F=l(.10498-5.37899*a),C=l(.42681-.41855*a),M=l(.14943-5.37511*a);var E=.84*z+.31*B+14.27*K+7.26*D+.28*F+.24*C;var J=2.94*z+.31*B+14.27*K+9.34*D+1.12*F+.83*C;var G=-6.4*z-1.89*C;K=.21*z+.31*B+14.27*K-88.7*D-15.3*F+.24*C-1.86*M;B=E-G;z=-3.332E-6*l(.59734-5.37261*a)-5.39E-7*l(.35498-5.37899*a)-6.4E-8*l(.39943-5.37511*a);E=H*t(.60643382+1336.85522467*a-3.13E-6*u)+E/L;J=H*t(.37489701+1325.55240982*a+2.565E-5*u)+J/L;G=H*t(.99312619+99.99735956*a-4.4E-7*u)+G/L;K=H*t(.25909118+1342.2278298*a-8.92E-6*u)+ -K/L;u=H*t(.82736186+1236.85308708*a-3.97E-6*u)+B/L;for(n=1;4>=n;++n){switch(n){case 1:var N=J;var x=4;var O=1.000002208;break;case 2:N=G;x=3;O=.997504612-.002495388*a;break;case 3:N=K;x=4;O=1.000002708+139.978*z;break;case 4:N=u,x=6,O=1}g(0,1,1);g(1,n,Math.cos(N)*O);h(0,n,0);h(1,n,Math.sin(N)*O);for(p=2;p<=x;++p)e(c(w,p-1,n),c(r,p-1,n),c(w,1,n),c(r,1,n),function(a,b){return g(p,n,a),h(p,n,b)});for(p=1;p<=x;++p)g(-p,n,c(w,p,n)),h(-p,n,-c(r,p,n))}f(13.902,14.06,-.001,.2607,0,0,0,4);f(.403,-4.01,.394, +K/L;u=H*t(.82736186+1236.85308708*a-3.97E-6*u)+B/L;for(n=1;4>=n;++n){switch(n){case 1:var N=J;var x=4;var O=1.000002208;break;case 2:N=G;x=3;O=.997504612-.002495388*a;break;case 3:N=K;x=4;O=1.000002708+139.978*z;break;case 4:N=u,x=6,O=1}g(0,n,1);g(1,n,Math.cos(N)*O);h(0,n,0);h(1,n,Math.sin(N)*O);for(p=2;p<=x;++p)e(c(w,p-1,n),c(r,p-1,n),c(w,1,n),c(r,1,n),function(a,b){return g(p,n,a),h(p,n,b)});for(p=1;p<=x;++p)g(-p,n,c(w,p,n)),h(-p,n,-c(r,p,n))}f(13.902,14.06,-.001,.2607,0,0,0,4);f(.403,-4.01,.394, .0023,0,0,0,3);f(2369.912,2373.36,.601,28.2333,0,0,0,2);f(-125.154,-112.79,-.725,-.9781,0,0,0,1);f(1.979,6.98,-.445,.0433,1,0,0,4);f(191.953,192.72,.029,3.0861,1,0,0,2);f(-8.466,-13.51,.455,-.1093,1,0,0,1);f(22639.5,22609.07,.079,186.5398,1,0,0,0);f(18.609,3.59,-.094,.0118,1,0,0,-1);f(-4586.465,-4578.13,-.077,34.3117,1,0,0,-2);f(3.215,5.44,.192,-.0386,1,0,0,-3);f(-38.428,-38.64,.001,.6008,1,0,0,-4);f(-.393,-1.43,-.092,.0086,1,0,0,-6);f(-.289,-1.59,.123,-.0053,0,1,0,4);f(-24.42,-25.1,.04,-.3,0,1,0, 2);f(18.023,17.93,.007,.1494,0,1,0,1);f(-668.146,-126.98,-1.302,-.3997,0,1,0,0);f(.56,.32,-.001,-.0037,0,1,0,-1);f(-165.145,-165.06,.054,1.9178,0,1,0,-2);f(-1.877,-6.46,-.416,.0339,0,1,0,-4);f(.213,1.02,-.074,.0054,2,0,0,4);f(14.387,14.78,-.017,.2833,2,0,0,2);f(-.586,-1.2,.054,-.01,2,0,0,1);f(769.016,767.96,.107,10.1657,2,0,0,0);f(1.75,2.01,-.018,.0155,2,0,0,-1);f(-211.656,-152.53,5.679,-.3039,2,0,0,-2);f(1.225,.91,-.03,-.0088,2,0,0,-3);f(-30.773,-34.07,-.308,.3722,2,0,0,-4);f(-.57,-1.4,-.074,.0109, 2,0,0,-6);f(-2.921,-11.75,.787,-.0484,1,1,0,2);f(1.267,1.52,-.022,.0164,1,1,0,1);f(-109.673,-115.18,.461,-.949,1,1,0,0);f(-205.962,-182.36,2.056,1.4437,1,1,0,-2);f(.233,.36,.012,-.0025,1,1,0,-3);f(-4.391,-9.66,-.471,.0673,1,1,0,-4);f(.283,1.53,-.111,.006,1,-1,0,4);f(14.577,31.7,-1.54,.2302,1,-1,0,2);f(147.687,138.76,.679,1.1528,1,-1,0,0);f(-1.089,.55,.021,0,1,-1,0,-1);f(28.475,23.59,-.443,-.2257,1,-1,0,-2);f(-.276,-.38,-.006,-.0036,1,-1,0,-3);f(.636,2.27,.146,-.0102,1,-1,0,-4);f(-.189,-1.68,.131, diff --git a/source/python/astronomy.py b/source/python/astronomy.py index 44ceb4c8..66df166a 100644 --- a/source/python/astronomy.py +++ b/source/python/astronomy.py @@ -675,3 +675,113 @@ def _geo_pos(time, observer): pos2 = _nutation(time, -1, pos1) outpos = _precession(time.tt, pos2, 0.0) return outpos + +def _spin(angle, pos1): + angr = angle * _DEG2RAD + cosang = math.cos(angr) + sinang = math.sin(angr) + return [ + +cosang*pos1[0] + sinang*pos1[1], + -sinang*pos1[0] + cosang*pos1[1], + pos1[2] + ] + +def _ter2cel(time, vec1): + gast = _sidereal_time(time) + return _spin(-15.0 * gast, vec1) + +#---------------------------------------------------------------------------- +# CalcMoon + +class _Array1: + def __init__(xmin, xmax): + self.min = xmin + self.array = [0] * (xmax - xmin + 1) + + def __getitem__(self, key): + return self.array[key - self.min] + + def __setitem__(self, key, value): + self.array[key - self.min] = value + +class _Array2: + def __init__(xmin, xmax, ymin, ymax): + self.min = xmin + self.array = [_Array1(ymin, ymax) for i in range(xmax - xmin + 1)] + + def __getitem__(self, key): + return self.array[key - self.min] + + def __setitem__(self, key, value): + self.array[key - self.min] = value + +def _CalcMoon(time): + T = time.tt / 36525 + co = _Array2(-6, 6, 1, 4) + si = _Array2(-6, 6, 1, 4) + + def AddThe(c1, s1, c2, s2): + return (c1*c2 - s1*s2, s1*c2 + c1*s2) + + def Sine(phi): + return math.sin(_PI2 * phi) + + def Frac(x): + return x - math.floor(x) + + T2 = T*T + DLAM = 0 + DS = 0 + GAM1C = 0 + SINPI = 3422.7000 + S1 = Sine(0.19833+0.05611*T) + S2 = Sine(0.27869+0.04508*T) + S3 = Sine(0.16827-0.36903*T) + S4 = Sine(0.34734-5.37261*T) + S5 = Sine(0.10498-5.37899*T) + S6 = Sine(0.42681-0.41855*T) + S7 = Sine(0.14943-5.37511*T) + DL0 = 0.84*S1+0.31*S2+14.27*S3+ 7.26*S4+ 0.28*S5+0.24*S6 + DL = 2.94*S1+0.31*S2+14.27*S3+ 9.34*S4+ 1.12*S5+0.83*S6 + DLS =-6.40*S1 -1.89*S6 + DF = 0.21*S1+0.31*S2+14.27*S3-88.70*S4-15.30*S5+0.24*S6-1.86*S7 + DD = DL0-DLS + DGAM = ((-3332E-9 * Sine(0.59734-5.37261*T) + -539E-9 * Sine(0.35498-5.37899*T) + -64E-9 * Sine(0.39943-5.37511*T))) + + L0 = PI2*Frac(0.60643382+1336.85522467*T-0.00000313*T2) + DL0/ARC + L = PI2*Frac(0.37489701+1325.55240982*T+0.00002565*T2) + DL /ARC + LS = PI2*Frac(0.99312619+ 99.99735956*T-0.00000044*T2) + DLS/ARC + F = PI2*Frac(0.25909118+1342.22782980*T-0.00000892*T2) + DF /ARC + D = PI2*Frac(0.82736186+1236.85308708*T-0.00000397*T2) + DD /ARC + + I = 1 + while I <= 4: + if I == 1: + ARG=L; MAX=4; FAC=1.000002208 + elif I == 2: + ARG=LS; MAX=3; FAC=0.997504612-0.002495388*T + elif I == 3: + ARG=F; MAX=4; FAC=1.000002708+139.978*DGAM + else: + ARG=D; MAX=6; FAC=1.0 + + co[0][I] = 1 + co[1][I] = math.cos(ARG) * FAC + si[0][I] = 0 + si[1][I] = math.sin(ARG) * FAC + + J = 2 + while J <= MAX: + co[J][I], si[J][I] = AddThe(CO(J-1,I), SI(J-1,I), CO(1,I), SI(1,I)) + J += 1 + + J = 1 + while J <= MAX: + co[-J][I] = +co[J][I] + si[-J][I] = -si[J][I] + J += 1 + + I += 1 +