github.com/pkimpel.retro-b5500
1
0
Fork 0
mirror of https://github.com/pkimpel/retro-b5500.git synced 2026-04-25 03:45:34 +00:00
Code Releases Activity
Files
d1e59a2f7c748e42aebee2aee60810fbe46be54a
pkimpel.retro-b5500/tests/E-mode/YUSPAR-XEM.alg_m

339 lines
30 KiB
Plaintext
Raw Blame History

% FREE STANDING YU/SPARROW ROUTINE -- 3/5/70 -- P H KIMPEL 00337000
BEGIN 00339000
FILE YUSPAR (KIND=READER, DEPENDENTSPECS); 00340000
FILE OUTPUT (KIND=PRINTER, MAXRECSIZE=132, FRAMESIZE=8, FILEUSE=OUT); 00341000
FORMAT FT ( "DATE: "A5" PROC TIME: "F8.2" I/O TIME: "F8.2" ELAPSED: " 00342000
F8.2" SECS"), 00343000
FMTYS (10(R20.12,/)); 00344000
INTEGER 00350000
STARTT, 00353000
NPOINTS, 00356000
NA, 00357000
NB, 00358000
NAB, 00359000
N1 , 00360000
N2, 00361000
N12, 00362000
KUSTAR, 00363000
NUSTAR, 00364000
NSUSTAR; 00365000
REAL 00366000
FN, 00367000
GMDT, 00368000
D, 00369000
H, 00370000
R, 00371000
MU1, 00372000
MU2, 00373000
DPDZ, 00374000
PI, 00375000
MU12, 00376000
ISUM, 00377000
YSTRT, 00378000
YSTEP, 00379000
YEND, 00380000
XSTEP, 00381000
XEND, 00382000
XX, 00383000
YY, 00384000
USUM, 00385000
Q1Q1FULL, 00386000
Q2Q2FULL, 00387000
Q1STARQ1FULL, 00387100
Q2STARQ2FULL, 00387200
Q1Q1STAR, 00387300
Q2Q2STAR, 00387400
X, 00387500
Y; 00387600
REAL TAB21, TAB22, Q1Q2; 00388000
00389000
% BEGIN EXECUTION IN OUTER BLOCK. 00390000
NPOINTS:= 100; NA:= NB:= 17; NAB:= NA + NB; 00391000
PI:= 3.141592653589796; 00403000
R:= 1.0; 00405000
FOR MU12:= 1.0,2.0,4.0,7.0,10.0,20.0,40.0,70.0,100.0,200.0,400.0, 00410000
700.0,1000.0,2000.0 DO 00411000
FOR H:= -0.8 STEP 0.2 UNTIL 0.8 DO 00412000
BEGIN 00413000
STARTT:= TIME(1); 00414000
%***********************************************************************00418000
BEGIN % ***** INNER BLOCK ***** 00419000
REAL ARRAY 00420000
W[1:NPOINTS], 00421000
Y[1:1, 1:NPOINTS+NAB], 00422000
Z[1:NAB, 1:NPOINTS+NAB], 00423000
X[1:NAB+1, 1:NPOINTS+NAB], 00423100
DEV[1:1, 1:NPOINTS], 00423200
COF[1:1, 1:NAB], 00423300
STD[1:1], 00423400
CV[1:NAB+1, 1:NAB], 00423500
VCV[1:1, 1:NAB+1, 1:NAB], 00423600
Q[1:1, 1:NAB+1], 00423700
Q2, 00423800
E, 00423900
EP[1:1, 1:NAB], 00424000
A[1:NAB, 1:NAB], 00424100
GF[1:NAB+1], 00424200
ENF[1:NAB]; 00424300
ARRAY COF1[1] = COF[1,*]; 00424400
LABEL DUMPIT; 00424500
DUMP OUTPUT (D, MU1, MU2, H, MU12, N1, N2, GMDT, Q1Q1STAR, Q2Q2STAR, 00424600
COF1, STD, TAB21, TAB22, Q1Q2) DUMPIT:1; 00424700
00425000
%***********************************************************************00426000
REAL PROCEDURE G (I, X, Y); 00427000
VALUE I, X, Y; REAL X, Y; INTEGER I; 00428000
CASE (I - 1) OF 00429000
BEGIN 00430000
G:= 1; % 100431000
G:= X; % 200432000
G:= Y; % 300433000
G:= X**2 - Y**2; % 400434000
G:= 2.0 * X * Y; % 500435000
G:= X**3 - 3.0 * X * Y**2; % 600436000
G:= 3.0 * X**2 * Y - Y**3; % 700437000
G:= X**4 + Y**4 - 6.0 * X**2 * Y**2; % 800438000
G:= 4.0 * X**3 * Y - 4.0 * X * Y**3; % 900439000
G:= X**5 - 10.0 * X**3 * Y**2 + 5.0 * X * Y**4; %1000440000
G:= Y**5 - 10.0 * Y**3 * X**2 + 5.0 * Y * X**4; %1100441000
G:= X**6 - 15.0 * X**4 * Y**2 + 15.0 * X**2 * Y**4 - Y**6; %1200442000
G:= 6.0 * X**5 * Y + 6.0 * X * Y**5 - 20.0 * X**3 * Y**3; %1300443000
G:= X**7 - 21.0 *X**5 * Y**2 + 35.0 * X**3 * Y**4 - 7.0 * X * %1400444000
Y**6; 00444100
G:= 7.0 *X**6 * Y - 35.0 * X**4 * Y**3 + 21.0 * X**2 * Y**5 - %1500445000
Y**7; 00445100
G:= X**8 + Y**8 - 28.0 *X**6 * Y**2 + 70.0 * X**4 * Y**4 - 28.0 00446000
* X**2 * Y**6; %1600447000
G:= 8.0 * X**7 *Y - 56.0 * X**5 * Y**3 + 56.0 * X**3 * Y**5 00448000
- 8.0 *X * Y**7; %1700449000
END G; 00450000
00451000
%***********************************************************************00452000
REAL PROCEDURE DG (I, X, Y); 00453000
VALUE I, X, Y; INTEGER I; REAL X, Y; 00454000
CASE (I - 1) OF 00455000
BEGIN 00456000
DG:= 0.0; % 100457000
DG:= 0.0; % 200458000
DG:= 1.0; % 300459000
DG:= -2.0 * Y; % 400460000
DG:= 2.0 * X; % 500461000
DG:= -6.0 * X *Y; % 600462000
DG:= 3.0 * X**2 - 3.0 * Y**2; % 700463000
DG:= 4.0 * Y**3 - 12.0 *X**2 * Y; % 800464000
DG:= 4.0 * X**3 - 12.0 * X * Y**2; % 900465000
DG:= -20.0 *X**3 * Y + 20.0 * X * Y**3; %1000466000
DG:= 5.0 * Y**4 - 30.0 * Y**2 *X**2 + 5.0 * X**4; %1100467000
DG:= -30.0 * X**4 *Y + 60.0 * X**2 * Y**3 - 6.0 * Y**5; %1200468000
DG:= 6.0 * X**5 + 30.0 * X * Y**4 - 60.0 * X**3 * Y**2; %1300469000
DG:= -42.0 * X**5 * Y + 140.0 * X**3 * Y**3 - 42.0 * X * Y**5;%1400470000
DG:= 7.0 * X**6 - 105.0 * X**4 * Y**2 + 105.0 * X**2 * Y**4 00471000
- 7.0 * Y**6; %1500472000
DG:= 8.0 *Y**7 - 56.0 * X**6 * Y + 280.0 * X**4 * Y**3 00473000
- 168.0 * X**2 * Y**5; %1600474000
DG:= 8.0 * X**7 - 168.0 * X**5 * Y**2 + 280.0 * X**3 * Y**4 00475000
- 56.0 * X * Y**6; %1700476000
END DG; 00477000
00478000
%***********************************************************************00479000
REAL PROCEDURE USTAR (I, X, Y); 00480000
VALUE I, X, Y; INTEGER I; REAL X, Y; 00481000
BEGIN 00482000
NUSTAR:= IF I = 1 THEN NA ELSE NA + NB; 00483000
NSUSTAR:= IF I = 1 THEN 1 ELSE NA + 1; 00484000
USUM:= 0.0; 00485000
FOR KUSTAR:= NSUSTAR STEP 1 UNTIL NUSTAR DO 00486000
USUM:= USUM + COF[1,KUSTAR] * G(KUSTAR-NSUSTAR+1, X, Y); 00487000
USTAR:= USUM; 00488000
END USTAR; 00489000
00490000
%***********************************************************************00491000
REAL PROCEDURE INTEGRALUSTAR (I, S); 00492000
VALUE I, S; INTEGER I; BOOLEAN S; 00493000
BEGIN 00494000
ISUM:= 0.0; 00495000
IF I = 1 THEN BEGIN YSTRT:= H; YEND:= 0.9999; END 00496000
ELSE BEGIN YSTRT:= -0.9999; YEND:= H; END; 00497000
YSTEP:= (YEND - YSTRT)/30.00; 00498000
FOR YY:= YSTRT STEP YSTEP UNTIL YEND DO 00499000
BEGIN XSTEP:= (XX:= SQRT(1.0 - YY**2))/15.0; 00500000
FOR XX:= -XX STEP XSTEP UNTIL 0.0 DO 00501000
ISUM:= ISUM + XSTEP * YSTEP * (USTAR (I, XX, YY) 00502000
- (IF S THEN 0.25 * (XX**2 + YY**2) ELSE 0.0)); 00503000
END; 00504000
INTEGRALUSTAR:= 2.0 * ISUM; 00505000
END INTEGRALUSTAR; 00506000
00507000
%***********************************************************************00508000
REAL PROCEDURE ARCSIN (X); 00509000
VALUE X; REAL X; 00510000
ARCSIN:= ARCTAN(X/SQRT(1.0 - X**2)); 00511000
00512000
%***********************************************************************00513000
REAL PROCEDURE EQUATION23 (A); 00514000
VALUE A; REAL A; 00515000
EQUATION23:= 0.5 - 2.0/3.0/PI * (A * (1.0 - A**2)**1.5 00516000
+ 1.5 * A *SQRT(1.0 - A**2) + 1.5 * ARCSIN(A)); 00517000
00518000
%***********************************************************************00519000
00520000
PROCEDURE ORTHO (W, Y, Z, N, M, P, COF); 00521000
VALUE N, M, P; 00522000
REAL ARRAY COF, W[1]; 00523000
REAL ARRAY Y, Z [1, 1]; 00524000
INTEGER N, M, P; 00525000
% 00526000
%COMMENT 00527000
% -ORTHO- IS TAKEN FROM ACM ALGORITHM 127 [COMM. ACM, VOL.5, 00528000
% OCTOBER 1962, P. 511, AUTHOR: PHILIP J. WALSH], 00529000
% AND INCORPORATES THE SUGGESTIONS IN BARRODALE-S 00530000
% CERTIFICATION [COMM. ACM, VOL.13, #2, FEB.1970, 122] 00531000
% .... P H KIMPEL 3/25/70 U OF DELAWARE 00534000
% 00535000
BEGIN 00536000
INTEGER NPP, M1, N2, M2, R1, RBAR, P2, BEI, I, J, K, THI, ALI, 00537000
OMI, R; 00538000
REAL ARRAY PK, XP [1:N+P], QK[1:M+1], X[1:M+1, 1:N+P]; 00539000
REAL SUM, DK2, DK; 00540000
LABEL BOX1, BOX2, BOX3, BOX4, BOX5, BOX6, BOX7, OM1, OM2, BOX8, BOX9, 00541000
TH1, TH2, TH3, BOX10, AL1, AL2, BOX12, BE1, BE2, BOX14, FINAL; 00542000
SWITCH BE:= BE1, BE2; 00543000
SWITCH TH:= TH1, TH2, TH3; SWITCH AL:= AL1, AL2; 00544000
SWITCH OM:= OM1, OM2; 00545000
NPP:= N + P; R:= 1; M1:= M-1; N2:= N+1; M2:= M+1; 00546000
R1:= 0; RBAR:= R; P2:= P+1; BEI:= 1; 00547000
BOX1: FOR I:= 1 STEP 1 UNTIL M DO BEGIN 00548000
FOR J:= 1 STEP 1 UNTIL N DO 00549000
X[I,J]:= Z[I,J] END; 00550000
BOX2: FOR I:= 1 STEP 1 UNTIL M DO BEGIN 00551000
FOR J:= N2 STEP 1 UNTIL NPP DO 00552000
X[I,J]:= 0.0; X[I,N+I]:= 1.0 END; GO TO BOX3; 00553000
BOX3: K:= 1; 00554000
BOX4: THI:= 1; 00555000
BOX5: ALI:= OMI:= 1; 00556000
FOR J:= 1 STEP 1 UNTIL P DO PK[N+J]:= 0.0; 00557000
BOX6: FOR I:= 1 STEP 1 UNTIL N DO 00558000
PK[I]:= X[K,I] * W[I]; 00559000
BOX7: GO TO OM[OMI]; 00560000
OM1: FOR I:= 1 STEP 1 UNTIL K DO BEGIN SUM:= 0.0; 00561000
FOR J:= 1 STEP 1 UNTIL NPP DO 00562000
SUM:= SUM + PK[J] * X[I,J]; QK[I]:= SUM END; 00563000
GO TO BOX8; 00564000
OM2: DK2:= 0.0; FOR I:= 1 STEP 1 UNTIL NPP DO 00565000
DK2:= DK2 + PK[I] * X[K,I]; 00566000
DK:= SQRT(DK2); 00567000
FOR I:= 1 STEP 1 UNTIL NPP DO 00568000
X[K,I]:= X[K,I]/DK;; 00569000
OMI:= 1; GO TO BOX6; 00570000
BOX8: FOR I:= 1 STEP 1 UNTIL K-1 DO 00571000
QK[I]:= -QK[I]; QK[K]:= 1.0; 00572000
FOR I:= 1 STEP 1 UNTIL NPP DO BEGIN 00573000
SUM:= 0.0; FOR J:= 1 STEP 1 UNTIL K DO 00574000
SUM:= SUM + X[J,I] * QK[J]; 00575000
XP[I]:= SUM END; GO TO BOX9; 00576000
BOX9: GO TO TH[THI]; 00577000
TH1: FOR I:= 1 STEP 1 UNTIL NPP DO 00578000
X[K,I]:= XP[I]; GO TO BOX10; 00579000
TH2: FOR I:= 1 STEP 1 UNTIL P DO 00580000
COF[I]:= -XP[N+I]; THI:= 3; GO TO TH1; 00581000
TH3: GO TO BOX14; 00582000
BOX10: GO TO AL[ALI]; 00583000
AL1: OMI:= ALI:= 2; GO TO BOX6; 00584000
AL2: IF K < M THEN BEGIN K:= K + 1; GO TO BOX4; END 00585000
ELSE GO TO BOX12; 00586000
BOX12: GO TO BE[BEI]; 00587000
BE1: BEI:= THI:= 2; K:= K + 1; GO TO BOX14; 00588000
BE2: GO TO BOX14; 00589000
BOX14: IF RBAR = 0 THEN GO TO FINAL ELSE RBAR:= RBAR - 1; 00590000
R1:= R1 + 1; THI:= 2; 00591000
FOR I:= 1 STEP 1 UNTIL N DO 00592000
X[M2,I]:= Y[R1,I]; 00593000
FOR I:= 1 STEP 1 UNTIL P DO 00594000
X[M2,N+I]:= 0.0; GO TO BOX5; 00595000
FINAL: END ORTHO ; 00596000
%***********************************************************************00597000
00598000
00598400
%***********************************************************************00598500
%***********************************************************************00599000
% BEGIN EXECUTION (INNER BLOCK) %00600000
%***********************************************************************00601000
%***********************************************************************00601500
BEGIN 00602000
COMMENT INITALIZE & SET UP DATA FOR -ORTHO-; 00604000
INTEGER I, J; 00605000
REAL THETA1, ARCC1, ARCC2, ARCC12, DS1, DS2, DS12, XX, YY, TOTALARC; 00606000
%% CALCULATE ARC LENGTHS OF EACH INTERFACE. 00607000
ARCC12:= 2.0 *SQRT(1.0 - H**2); 00608000
ARCC1:= 2.0 * (THETA1:= PI/2.0 - ARCTAN(2.0*H/ARCC12)); 00609000
ARCC2:= 2.0 *PI - ARCC1; 00610000
%% CALCULATE NO. OF POINTS ON EACH ARC. 00611000
N1:= NPOINTS * ARCC1/(TOTALARC:= ARCC1 + ARCC2 + 2*ARCC12); 00612000
N2:= NPOINTS *ARCC2/TOTALARC; 00613000
IF (N12:= NPOINTS - N1 - N2)MOD 2 = 0 THEN 00614000
N12:= N12 DIV 2 00615000
ELSE BEGIN 00616000
N12:= N12 DIV 2 + 1; 00617000
N2:= N2 - 1; 00618000
END; 00619000
%% CALCULATE THE LENGTH OF LINE SEGMENTS BETWEEN POINTS ON EACH 00620000
%% INTERFACE. 00621000
DS1:= ARCC1/(N1 - 1); 00622000
DS2:= ARCC2/(N2 - 1); 00623000
DS12:= ARCC12/(N12 - 1); 00624000
00625000
FOR J:= 1 STEP 1 UNTIL NPOINTS-N12 DO 00626000
IF J <= N1 THEN %% DATA FOR ARC C1. 00627000
BEGIN 00628000
W[J]:= 1.0; Y[1,J]:= 0.25; 00629000
XX:= SIN(YY:= DS1*(J-1)-THETA1); YY:= COS(YY); 00630000
FOR I:= 1 STEP 1 UNTIL NAB DO 00631000
Z[I,J]:= IF I <= NA THEN G(I, XX, YY) ELSE 0.0; 00632000
END ELSE 00633000
IF J <= N1 + N2 THEN %% DATA FOR ARC C2. 00634000
BEGIN 00635000
W[J]:= DS2/DS1; Y[1,J]:= 0.25; 00636000
XX:= SIN(YY:= - THETA1 - DS2*(J-N1-1)); YY:= COS(YY); 00637000
FOR I:= 1 STEP 1 UNTIL NAB DO 00638000
Z[I,J]:= IF I <= NA THEN 0.0 ELSE G(I-NA, XX, YY); 00639000
END ELSE %% DATA FOR ARC C12. 00640000
BEGIN 00641000
W[J]:= W[J+N12]:= DS12/DS1; 00642000
Y[1,J]:= 0.25*((XX:= DS12*(J-N1-N2-1)-ARCC12/2)**2 + H**2)00643000
* (1.0 - MU12); 00644000
Y[1,J+N12]:= 0.0; 00645000
FOR I:= 1 STEP 1 UNTIL NAB DO 00646000
IF I <= NA THEN 00647000
BEGIN 00648000
Z[I,J]:= G(I, XX, H); 00649000
Z[I,J+N12]:= DG(I, XX, H); 00650000
END ELSE 00651000
BEGIN 00652000
Z[I,J]:= -G (I-NA, XX, H) * MU12; 00653000
Z[I,J+N12]:= -DG (I-NA, XX, H); 00654000
END; 00655000
END; 00656000
END INITIALIZATION; 00657000
00658000
%%% THE NUMBER SMASHING COMMENCES. 00659000
ORTHO (W, Y, Z, NPOINTS, NAB, NAB, COF[1,*]); %######################00660000
00661000
%%% UPON EXIT FROM ORTHO, THE YU/SPARROW COEFFICIENTS, A[J], B[J], 00662000
%%%% ARE CONTAINED IN THE VECTOR COF[1,*]. 00663000
00664000
Q1Q1FULL:= 8.0/PI * INTEGRALUSTAR(1,TRUE); 00665000
Q2Q2FULL:= 8.0/PI * INTEGRALUSTAR(2,TRUE); 00666000
Q1STARQ1FULL:= EQUATION23(H); 00667000
Q2STARQ2FULL:= 1.0 - Q1STARQ1FULL; 00668000
TAB21:= Q1Q1FULL * PI/8.0; TAB22:= Q2Q2FULL * PI/8.0; 00669000
Q1Q2:= Q1Q1FULL/Q2Q2FULL/MU12; 00670000
Q1Q1STAR:= Q1Q1FULL/Q1STARQ1FULL; 00671000
Q2Q2STAR:= Q2Q2FULL/Q2STARQ2FULL; 00671500
00672000
00673000
DUMPIT: 00680000
WRITE(OUTPUT, FT, TIME(0), TIME(2)/60, TIME(3)/60, (TIME(1)-STARTT)/60);00681000
END INNER BLOCK; 00682000
END H FOR LOOP; 00684000
END. 00686000
View Git Blame Copy Permalink