|
DataMuseum.dkPresents historical artifacts from the history of: RC4000/8000/9000 |
This is an automatic "excavation" of a thematic subset of
See our Wiki for more about RC4000/8000/9000 Excavated with: AutoArchaeologist - Free & Open Source Software. |
top - metrics - download
Length: 7680 (0x1e00)
Types: TextFile
Names: »bshlmproctx«
└─⟦621cfb9a2⟧ Bits:30002817 RC8000 Dump tape fra HCØ. Detaljer om "HC8000" projekt.
└─⟦0364f57e3⟧
└─⟦80900d603⟧ »giprocfile«
└─⟦00964e8f7⟧ Bits:30007478 RC8000 Dump tape fra HCØ.
└─⟦b2ec5d50f⟧
└─⟦80900d603⟧ »giprocfile«
└─⟦this⟧
; bshlmproctx * page 1 26 09 77, 9.41;
; bshlm_proc
; **********
if listing.yes
char 10 12 10
bshlm_proc = set 1
bshlm_proc = algol
external long procedure bshlm_proc
__________________________________
_ (reg_lab, mode, lat, az_dlat, s_dlng, dlng_s, az_2);
value mode, lat, az_dlat, s_dlng;
array reg_lab;
integer mode;
long lat, az_dlat, s_dlng, dlng_s, az_2;
comment
bshlm_proc (return) long
result depending on mode.
mode = 0 : the difference of latitude in the solution of
the first geodetic main problem
mode > 0 : the outgoing azimuth from station 1 towards
station 2
reg_lab (call) array
the values of the semimajor axis of the ellipsoid as a
long in units of '-6 m in reg_lab(6), and the flattening as
a real in reg_lab(7). the array is not changed by the
procedure and no other elements are used
mode (call) integer
mode < 0 is interpreted as mode = 0.
mode = 0 specifies a solution of the first (or direct) main
problem.
mode > 0 specifies a solution of the second (or inverse)
main problem. mode=1 is the shortest line, mode = 2 is the
line the other way around the earth. all following odd
modes is mode//2 wrap-arounds of mode 1 and similarly for
the even modes and mode 2
lat (call) long
mode = 0 : the latitude of the station 1
mode > 0 : the m_e_a_n latitude of the stations
az_dlat (call) long
mode = 0 : the outgoing azimuth
mode > 0 : the difference of latitude stat2 minus stat1
s_dlng (call) long
mode = 0 : the distance from station 1 to station 2
mode > 0 : the difference of longitude stat2 minus stat1
dlng_s (return) long
mode = 0 : the difference of longitude stat2 minus stat1
mode > 0 : the distance between the stat
az_2 (return) long
the azimuth at stat2 towards stat1;
\f
comment bshlmproctx * page 2 26 09 77, 9.41
0 1 2 3 4 5 6 7 8 9 ;
comment
all angles are geotype variables i.e. the interval -pi to pi
is mapped on the integer interval -2**47 to 2**47 - 1.
all distances are in units of '-6 m;
begin real b0, f, e2, n1, B1, A1, S, dB, dL, A2;
real procedure RED(a); value a; real a;
begin long w;
w:= a/rg;
RED := w*rg;
end RED;
real procedure BSHLM1(B1, A1, S, dB, dL, A2);
value B1, A1, S; real B1, A1, S, dB, dL, A2;
begin real bm, s, S1, SP, PF, k, b1r, b2r, B2, dl, a1, a2, s1sf;
array KES, KSE(1:3); integer i;
b1r:=red(B1);
bm:=bmax(b1r, A1, s1sf, PF, k, KSE);
SP:=S/PF;
KES(1):=-1+9*(k**2)/16; KES(2):=5*k/8; KES(3):=-29*(k**2)/48;
S1:=s1sf+0.5*k*sum(i, KSE(i)*sin(2*i*s1sf));
s:=SP+0.5*k*sum(i, KES(i)*cos(2*i*(S1+SP))*sin(2*i*SP));
B2:=0.5*pi-red(2.0*sftr(A1, 0.5*pi-b1r, 2.0*s, A2, dl));
dL:=RED( dl-2.0*lmbd(s1sf, s, A1, dl, bm, k));
BSHLM1:=dB:=B2-B1;
end BSHLM1;
real procedure BSHLM2(mode, B, dB, dL, A1, S, A2);
value mode, B, dB, dL; real B, dB, dL, S, A1, A2;
integer mode;
begin
real bm, b1r, b2r, s1sf, s, SF, PF,
dl, dl1, k, mm2pi, pimode, last_bm;
array KSE(1:3);
integer i, iter;
iter := 0;
pimode := (mode//2)*pi;
mode := mode - 1;
mm2pi := pi*(mode mod 2);
b1r:=red(B-0.5*dB); b2r:=red(B+0.5*dB); dL:=RED(0.5*dL);
for dl1:=dL, dL+lmbd(s1sf, s, A1, dl, bm, k)
while abs(dl - dl1) > (abs dl+ abs s)*'-10 and iter<=10 000,
dL + lmbd(s1sf, s, A1, dl, bm, k) do
begin
iter := iter + 1;
if iter mod 1000 = 0 then
begin
write(out, nl, 1);
writegeocr(out, dl1, false add 70);
write(out, sp, 4);
writegeocr(out, dl-dl1, false add 70);
end;
dl := if iter<10 then dl1
_ else (0.5*(dl + dl1));
s:=sftr(RED(2*dl), 0.5*pi-b1r, 0.5*pi-b2r, A2, A1);
s:=pimode + (if mode mod 2 = 0 then s else -s);
A1 := RED(A1 + mm2pi);
A2 := RED(A2 + mm2pi);
bm:=bmax(b1r, A1, s1sf, PF, k, KSE);
last_bm := bm := if iter<10 then bm
_ else (0.5*(bm + last_bm));
\f
comment bshlmproctx * page 3 26 09 77, 9.41
0 1 2 3 4 5 6 7 8 9 ;
if iter mod 1000=0 then
begin
write(out, nl, 1);
writegeocr(out, bm, false add 70);
end;
end;
S:=PF*(s+0.5*k*
sum(i, KSE(i)*cos(2.0*i*(s1sf+s))*sin(2.0*i*s)));
if iter>=10 000 then S := A1 := A2 := 0;
BSHLM2 := A1;
end BSHLM2;
real procedure sftr(C, a, b, A, B);
value C, a, b; real A, B, C, a, b;
begin real svscs, svdss, cvscs, cvdss, cosCh, sinCh, sum, dif;
sum:=0.5*(a+b); dif:=0.5*(a-b); C:=0.5*C;
cosCh:=cos(C); sinCh:=sin(C);
svscs:=cosCh*cos(dif); svdss:=cosCh*sin(dif);
cvscs:=sinCh*cos(sum); cvdss:=sinCh*sin(sum);
sum:=arg(cvscs, svscs); dif:=arg(cvdss, svdss);
A:=RED(-sum-dif); B:=RED(sum-dif);
sftr:=RED(
arg(sqrt(svscs**2+cvscs**2), sqrt(svdss**2+cvdss**2)));
end sftr;
real procedure red(B);
value B; real B;
red:=arg(cos(B), (1-f)*sin(B));
real procedure bmax(b1, A1, s1sf, PF, k, KSE);
value b1, A1; real b1, A1, s1sf, PF, k; array KSE;
begin real c1, sc, ss, cosb1, bm, k1;
cosb1:=cos(b1); c1:=abs(sin(A1)*cosb1);
sc:=sin(b1); ss:=-cosb1*cos(A1);
s1sf:=arg(sc, ss);
bmax:=bm:=arg(c1, sqrt(ss**2+sc**2));
k1:=e2*sin(bm)**2/(1-e2*cos(bm)**2)/4.0;
k:=k1:=k1*(1+k1*(2.0+k1*(5.0+k1*14.0)));
PF:=b0*(2.0+0.5*k1**2)/(1-k1);
KSE(1):=1-3*(k1**2)/8; KSE(2):=-k1/8; KSE(3):=(k1**2)/24;
end bmax;
real procedure sum(i, x); integer i; real x;
begin real s;
s:=0; for i:=3 step -1 until 1 do s:=s+x;
sum:=s;
end sum;
real procedure lmbd(s1sf, s, A1, dl, bm, k);
value s1sf, s, A1, dl, bm, k; real s1sf, s, A1, dl, bm, k;
begin real ss;
ss:=2*(s1sf + s);
lmbd := (if A1 >= 0 then 0.5 else -0.5)
*e2*cos(bm)*((n1-k/2*(1+k/2))*s
-k/4*(cos(ss)*sin(2*s)-k/4*cos(2*ss)*sin(4*s)));
end lmbd;
\f
comment bshlmproctx * page 4 26 09 77, 9.41
0 1 2 3 4 5 6 7 8 9 ;
f := reg_lab(7);
b0 := (long reg_lab(6))*'-6*(1-f);
e2 := f*(2-f);
n1 := 1 + f/(2-f);
B1 := lat*rg;
if mode <= 0 then
begin comment direct problem;
A1 := az_dlat*rg;
S := s_dlng*'-6;
bshlm_proc := BSHLM1(B1, A1, S, dB, dL, A2)/rg;
dlng_s := dL/rg;
az_2 := A2/rg;
end direct probelm
else
begin comment inverse problem;
dB := az_dlat*rg;
dL := s_dlng*rg;
bshlm_proc := BSHLM2(mode, B1, dB, dL, A1, S, A2)/rg;
dlng_s := S*'6;
az_2 := A2/rg;
end inverse problem;
end bshlm_proc;
end
if warning.yes
(mode 0.yes
message bshlm_proc not ok
lookup bshlmproc)
▶EOF◀