-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathtstAstro.en
143 lines (143 loc) · 6.08 KB
/
tstAstro.en
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
% * Grundlagen
>lat:= 47.5°; load "astro.e";
>now=day(2022,1,13,12,21,30)
8048.01458333
>printday(now), getymdhms(now)
2022-01-13 12:21:00
[2022, 1, 13, 12, 21, 0]
>locMünchenstein
[7.61322, 47.5212, 335.78]
>here=[locMünchenstein,5,1013];
>sun(now)
[295.129, -21.4295, 0.983523]
>raltaz(now,here,sun(now))
[190.624, 20.4255]
% * Sonnenauf- und untergänge
>function SunAz(t=now,loc=here) ....
$h=raltaz(t,loc,sun(t));
$return h[1];
$endfunction
>function SunAlt(t=now,loc=here)...
$h=raltaz(t,loc,sun(t));
$$return h[2];
$endfunction
>SunAz(), SunAz(day(2022,3,25,11,35.5))
190.624207618
180.000996406
>printday(rise("sun",now,here)), printday(set("sun",now,here))
2022-01-14 07:15:08
2022-01-13 16:01:00
>plot2d("SunAzmap(daymap(2022,6,21,x))",a=3.5,b=19.5,color=2, grid=3);
>plot2d("SunAzmap(daymap(2022,3,21,x))",color=3,add=1);
>plot2d("SunAzmap(daymap(2022,12,21,x))",color=12,add=1);
>plot2d("SunAzmap(daymap(2022,9,23,x))",color=4,add=1);
>xgrid(0:1:24); ygrid([60:30:300]); ygrid([103,283],color=1); ygrid([54,126,234,306],color=2); insimg(40);
%image% tstAstro-001.png
>plot2d("SunAltmap(daymap(2022,6,21,x))",a=3.5,b=19.5,c=-15,d=75,color=2,grid=1);
>plot2d("SunAltmap(daymap(2022,3,21,x))",color=3,add=1);
>plot2d("SunAltmap(daymap(2022,9,23,x))",color=4,add=1);
>plot2d("SunAltmap(daymap(2022,12,21,x))",color=12,add=1);
>xgrid(0:1:24); ygrid([-10:10:60]); ygrid([19,42.5,66],color=2); insimg(40);
%image% tstAstro-002.png
% Auf- und Untergänge über Zeitintervall (zZ deaktiviert!)
>..s=day(2022,1,1):1:day(2022,12,31); l=ones([2,cols(s)]); for i=1 to cols(s); l[1,i]:=rise("sun",s[i],here); end;
>..for i=1 to cols(s); l[2,i]:=set("sun",s[i],here); end;
>..raltaz(l[1,-1],here,sun(l[1,-1])), raltaz(l[2,-1],here,sun(l[2,-1]))
>..rA=(l[1]_SunAzmap(l[1])_l[2]_SunAzmap(l[2]))';
>..writematrix(rA,"riseAzi.txt");
>{MT,hd}=readtable("riseAzi.txt"); writetable(MT[1:8],dc=5,wc=15,fixed=[1,0,1,0],labc=["tRise","AziRise","tSet","AziSet"]);
tRise AziRise tSet AziSet
8036.80496 124.418 8037.15840 235.635
8037.80491 124.262 8038.15911 235.796
8038.80483 124.094 8039.15984 235.968
8039.80471 123.914 8040.16059 236.153
8040.80456 123.722 8041.16137 236.349
8041.80437 123.518 8042.16217 236.558
8042.80416 123.302 8043.16298 236.778
8043.80391 123.075 8044.16382 237.01
% Erstellen von Planetentabellen
>function SunPos(t1,t2,inc)...
$global here
$c=1;
$suntab=zeros([(t2-t1)/inc,3]);
$for t=t1 to t2 step inc;
$pos=raltaz(t,here,sun(t));
$suntab(c,1)=t; suntab(c,2)=pos(1); suntab(c,3)=pos(2);
$c=c+1;
$end;
$return suntab
$endfunction
% * Vertikale Sonnenuhr
% Berechnungselemente, Besonnungsdauer und Streiflicht
>t=secant("SunAz(day(2022,4,20,x))-283",17,19); hmsprint(t), raltaz(day(2022,4,20,t),here,sun(day(2022,4,20,t)))
17:56:14
[283, 4.16796]
>t=secant("SunAz(day(2022,4,20,x))-103",6,9); hmsprint(t), raltaz(day(2022,4,20,t),here,sun(day(2022,4,20,t)))
07:19:07
[103, 27.0808]
% Streiflichtberechnung
>d=day(2022,2,22,6):1:day(2022,10,19,6);
>l=ones([2,cols(d)]); for i=1 to cols(d); t=secant("SunAz(x,here)-103",d[i]); l[1,i]=t; l[2,i]=raltaz(t,here,sun(t))[2]; end;
>writematrix(l',"Streiflicht.txt");
>{MT,hd}=readtable("Streiflicht.txt"); writetable(MT[1:5],dc=5,wc=15,fixed=[1,0],labc=["day","Höhe"]);
day Höhe
8088.76221 -0.7271
8089.76305 -0.30815
8090.76388 0.10325
8091.76470 0.52445
8092.76552 0.9604
>printday(MT[1,1]), printday(MT[-1,1])
2022-02-23 06:17:35
2022-10-19 05:48:41
% geogr. Breite, Länge und Wandabweichung (West)
>phi=47.5212; la=-7.61322; wa=13.55;
% Erhebungs- und Substilarwinkel, Hilfsstundenwinkel
>g=dasin(dcos(phi)*dcos(wa)), f=datan(dsin(wa)/dtan(phi)), t0=datan(dsin(phi)/dtan(wa))
41.035040359
12.1082102207
71.9041517919
% Deklination bei Sonnenauf bzw. -untergang in Wandflucht
>del0:=asin(dsin(wa)*dcos(phi)); degprint(del0)
9°6'13.50''
% Komponenten des Schattenstabs
>{x,y,z}=rect(rad(f),rad(g),1); format(8,3); [0.2,0.25,0.5,1]'*[x,y,z], longformat;
0.148 0.032 0.131
0.184 0.040 0.164
0.369 0.079 0.328
0.738 0.158 0.657
>function bd(dekli, wa=wa) ...
$global phi;
$if wa==0 then break else
$t0=datan(dsin(phi)/dtan(wa));
$g=dasin(dcos(phi)*dcos(wa));
$return mod(t0+dasin(dtan(g)*dtan(dekli)),360);
$endif;
$endfunction
>tmp:=(bd(-23.44)_bd(0)_bd(23.44))'
[49.7338718293, 71.9041517919, 94.0744317545]
>bd(-23.44)
49.7338718293
>plot2d("bdmap(x,[10:10:50]')", -23.44, 23.44);
>plot2d("bd",add=1,color=2,thickness=2); xgrid([-35:5:35]); ygrid([0:15:90]); insimg(30);
%image% tstAstro-003.png
% * Analemma (Deklination / Zeitgleichung)
>t=day(2022,1,1):1:day(2022,12,31); ...
>l=ones([2,cols(t)]); for i=1 to cols(t); hi:=sun(t[i]); l[1,i]:=hi[2]; l[2,i]:=(-hi[1]+gst(t[i])); end; ...
>l[2]:=(mod(l[2],360)-180)/15*60; writematrix((t_l)',"Analemma.txt");
>aspect(0.5); plot2d(l[2],l[1],a=-15,b=17,c=-25,d=25,xl="Zeitgleichung [min]", yl="Dekl. [°]",>points,style=".",grid=0,>vertical); ...
>frame; xgrid(-15:5:15); ygrid(-25:5:25);...
>h:=ones([1,12]); for m=1 to 12; h[m]:=day(2016,m,21); end; ...
>h[9]:=h[9]+2; ...
>L=ones([2,cols(h)]); for i=1 to cols(h); hi:=sun(h[i]); L[1,i]:=hi[2]; L[2,i]:=(-hi[1]+gst(h[i])); end; ...
>L[2]:=(mod(L[2],360)-180)/15*60;
>plot2d(L[2],L[1], points=1, style="<>", add=1);
>h:=ones([1,12]); for m=1 to 12; h[m]:=day(2016,m,1); end; ...
>L=ones([2,cols(h)]); for i=1 to cols(h); hi:=sun(h[i]); L[1,i]:=hi[2]; L[2,i]:=(-hi[1]+gst(h[i])); end; ...
>L[2]:=(mod(L[2],360)-180)/15*60;
>plot2d(L[2],L[1], points=1, style="*", add=1);
>h:=ones([1,12]); for m=1 to 12; h[m]:=day(2016,m,11); end; ...
>L=ones([2,cols(h)]); for i=1 to cols(h); hi:=sun(h[i]); L[1,i]:=hi[2]; L[2,i]:=(-hi[1]+gst(h[i])); end; ...
>L[2]:=(mod(L[2],360)-180)/15*60;
>plot2d(L[2],L[1], points=1, style="o", add=1,); insimg; aspect();
%image% tstAstro-004.png
>