1
0
mirror of https://github.com/pkimpel/retro-b5500.git synced 2026-02-15 04:26:29 +00:00

Commit additional tests for B5500/E-mode floating point comparison tests exchanged with Sid Harg.

This commit is contained in:
paul.kimpel@digm.com
2012-08-29 16:53:06 +00:00
parent 691bc6be6e
commit 693998cb59
4 changed files with 1274 additions and 0 deletions

335
tests/ALGOLXEM/YUSPAR.alg_m Normal file
View File

@@ -0,0 +1,335 @@
% FREE STANDING YU/SPARROW ROUTINE -- 3/5/70 -- P H KIMPEL 00337000
BEGIN 00339000
FILE IN YUSPAR (1,10); 00340000
FILE OUT OUTPUT 18 "YU/SPAR" "RESULTS" (2,17); 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
LABEL DUMPIT; 00424500
DUMP OUTPUT (D, MU1, MU2, H, MU12, N1, N2, GMDT, Q1Q1STAR, Q2Q2STAR, 00424600
COF, 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 | Y*6; %1400444000
G ~ 7.0 |X*6 | Y - 35.0 | X*4 | Y*3 + 21.0 | X*2 | Y*5 - Y*7; %1500445000
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 File

@@ -0,0 +1,257 @@
% -ORTHO- TEST, 3/3/70 00000100120429PK
BEGIN 00000200
INTEGER STARTT; 00000300120429PK
FORMAT FT (//"DATE: "A6" PROC TIME: "F8.2" I/O TIME: "F8.2" ELAPSED: " 00000400120429PK
F8.2" SEC"); 00000500120429PK
INTEGER I, J, N, M, P, R, AI, AUI, ZEI, MUI; 00000600120429PK
REAL SUM, GMDT, FN; 00000700120429PK
FILE PARAMS (KIND=READER, DEPENDENTSPECS); 00000800120826PK
FILE FLUSH (KIND=PRINTER, MAXRECSIZE=132, FRAMESIZE=8, FILEUSE=OUT); 00000900120826PK
00001000120429PK
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%00002000120429PK
PROCEDURE ORTHO (W, Y, Z, N, FN, M, P, R, AI, AUI, MUI, ZEI, X, DEV, 00003000
COF, STD, CV, VCV, GMDT, Q, Q2, E, EP, A, GF, ENF); 00004000
VALUE N, M, P, R, AI, AUI, MUI, ZEI; 00005000
REAL FN, GMDT; 00006000
REAL ARRAY STD, GF, W, ENF [1]; 00007000
REAL ARRAY Y, Z, X, DEV, COF, CV, Q, Q2, E, EP, A [1, 1]; 00008000
REAL ARRAY VCV [1, 1, 1]; 00009000
INTEGER N, M, P, R, AI, AUI, ZEI, MUI; 00010000
00011000120429PK
COMMENT 00012000
ORTHO IS TAKEN FROM ACM ALGORITHM 127 [COMM. ACM, VOL.5, 00013000
OCTOBER 1962, P. 511, AUTHOR: PHILIP J. WALSH]; 00014000
00015000120429PK
BEGIN 00016000
INTEGER NPP, NPM, M1, N2, M2, R1, RBAR, P2, BEI, RHI, I18, GAI, SII, I,00017000
J, DEI, NUI, E1Z2, E1Z1, K, THI, ALI, OMI, NII; 00018000
REAL ARRAY PK, XP [1:N+P], QK[1:M+1]; 00019000
REAL DENOM, SUM, DK2, DK, FI, SS, SSQ; 00020000
LABEL BOX1, AT1, AT2, BOX2, AU1, AU2, BOX3, BOX4, BOX5, BOX6, MU1, 00021000
MU2, BOX7, OM1, OM2, BOX8, DE1, BOX8A, BOX8B, DE2, BOX8C, BOX8D, 00022000
BOX9, TH1, TH2, TH3, BOX10, AL1, AL2, BOX11, NU1, NU2, BOX12, 00023000
BE1, BE2, BOX13, GA1, GA2, BOX14, RH1, ZE1, ZE2, RH2, SI1, SI2, 00024000
FINAL; 00025000
SWITCH AT:= AT1, AT2; SWITCH ZE:= ZE1, ZE2; 00026000120826PK
SWITCH AU:= AU1, AU2; SWITCH MU:= MU1, MU2; 00027000120826PK
SWITCH BE:= BE1, BE2; SWITCH RH:= RH1, RH2; SWITCH GA:= GA1, GA2; 00028000120826PK
SWITCH SI:= SI1, SI2; SWITCH DE:= DE1, DE2; SWITCH NU:= NU1, NU2; 00029000120826PK
SWITCH TH:= TH1, TH2, TH3; SWITCH AL:= AL1, AL2; 00030000120826PK
SWITCH OM:= OM1, OM2; 00031000120826PK
NPP:= N+P; NPM:= N+M; M1:= M-1; N2:= N+1; M2:= M+1; 00032000120826PK
R1:= 0; RBAR:= R; P2:= P+1; DENOM:= IF N=M THEN 1.0 00033000120826PK
ELSE SQRT(N-M); BEI:= RHI:= I18:= 1; 00034000120826PK
IF (P ^= 0) THEN GAI:= SII:= 2 ELSE GAI:= SII:= 1; 00035000120826PK
BOX1: GO TO AT[AI]; 00036000
AT1: FOR J:= 1 STEP 1 UNTIL N DO BEGIN 00037000120826PK
X[2,J]:= Z[1,J]; X[1,J]:= 1.0 END; 00038000120826PK
FOR I:= 2 STEP 1 UNTIL M1 DO BEGIN 00039000120826PK
FOR J:= 1 STEP 1 UNTIL N DO 00040000120826PK
X[I+1,J]:= X[I,J] * X[2,J] END; GO TO BOX2; 00041000120826PK
AT2: FOR I:= 1 STEP 1 UNTIL M DO BEGIN 00042000120826PK
FOR J:= 1 STEP 1 UNTIL N DO 00043000120826PK
X[I,J]:= Z[I,J] END; 00044000120826PK
BOX2: IF P = 0 THEN GO TO BOX3 ELSE GO TO AU[AUI]; 00045000
AU1: FOR I:= 1 STEP 1 UNTIL M DO BEGIN 00046000120826PK
FOR J:= N2 STEP 1 UNTIL NPP DO 00047000120826PK
X[I,J]:= 0.0; X[I,N+I]:= 1.0 END; GO TO BOX3; 00048000120826PK
AU2: FOR I:= 1 STEP 1 UNTIL M DO BEGIN 00049000120826PK
FOR J:= N2 STEP 1 UNTIL NPP DO 00050000120826PK
X[I,J]:= Z[I,J] END; 00051000120826PK
BOX3: DEI:= NUI:= E1Z1:= E1Z2:= K:= 1; 00052000120826PK
BOX4: THI:= 1; 00053000120826PK
BOX5: ALI:= OMI:= 1; IF P = 0 THEN GO TO BOX6 ELSE 00054000120826PK
FOR J:= 1 STEP 1 UNTIL P DO PK[N+J]:= 0.0; 00055000120826PK
BOX6: GO TO MU[MUI]; 00056000
MU1: FOR I:= 1 STEP 1 UNTIL N DO PK[I]:= X[K,I]; 00057000120826PK
GO TO BOX7; 00058000
MU2: FOR I:= 1 STEP 1 UNTIL N DO 00059000120826PK
PK[I]:= X[K,I] * W[I]; GO TO BOX7; 00060000120826PK
BOX7: GO TO OM[OMI]; 00061000
OM1: FOR I:= 1 STEP 1 UNTIL K DO BEGIN SUM:= 0.0; 00062000120826PK
FOR J:= 1 STEP 1 UNTIL NPP DO 00063000120826PK
SUM:= SUM + PK[J] * X[I,J]; QK[I]:= SUM END; 00064000120826PK
GO TO BOX8; 00065000
OM2: DK2:= 0.0; FOR I:= 1 STEP 1 UNTIL NPP DO 00066000120826PK
DK2:= DK2 + PK[I] * X[K,I]; 00067000120826PK
DK:= SQRT(DK2); 00068000120826PK
GF[I18]:= DK; I18:= I18 + 1; 00069000120826PK
FOR I:= 1 STEP 1 UNTIL NPP DO 00070000120826PK
X[K,I]:= X[K,I]/DK;; 00071000120826PK
OMI:= 1; GO TO BOX6; 00072000120826PK
BOX8: GO TO DE[DEI]; 00073000
DE1: E1Z1:= -E1Z1; IF E1Z1 < 0 THEN GO TO BOX8B ELSE 00074000120826PK
GO TO BOX8A; 00075000
BOX8A: FOR I:= 1 STEP 1 UNTIL K-1 DO 00076000120826PK
QK[I]:= -QK[I]; QK[K]:= 1.0; 00077000120826PK
FOR I:= 1 STEP 1 UNTIL NPP DO BEGIN 00078000120826PK
SUM:= 0.0; FOR J:= 1 STEP 1 UNTIL K DO 00079000120826PK
SUM:= SUM + X[J,I] * QK[J]; 00080000120826PK
XP[I]:= SUM END; GO TO BOX9; 00081000120826PK
BOX8B: ENF[I18]:= SQRT(QK[K]); GO TO BOX8A; 00082000120826PK
DE2: E1Z2:= -E1Z2; IF E1Z2 < 0 THEN GO TO BOX8C ELSE 00083000120826PK
GO TO BOX8A; 00084000
BOX8C: FOR I:= 1 STEP 1 UNTIL M DO BEGIN 00085000120826PK
Q[R1,I]:= QK[I]; Q2[R1,I]:= QK[I] * QK[I] END; 00086000120826PK
Q[R1,M2]:= QK[M2]; E[R1,1]:= Q[R1,M2] - Q2[R1,1]; 00087000120826PK
FOR J:= 2 STEP 1 UNTIL M DO 00088000120826PK
E[R1,J]:= E[R1,J-1] - Q2[R1,J]; 00089000120826PK
FI:= 1.0; 00090000120826PK
FOR I:= 1 STEP 1 UNTIL M DO BEGIN 00091000120826PK
IF (FN - FI) > 0.0 THEN BEGIN IF E[R1,I] < 0.0 THEN 00093000
BEGIN EP[R1,I]:= -SQRT(ABS(E[R1,I])/(FN - FI)); 00094000120826PK
GO TO BOX8D; END 00095000
ELSE EP[R1,I]:= SQRT(E[R1,I]/(FN - FI)); 00096000120826PK
GO TO BOX8D; END ELSE E[R1,I]:= -1.0; 00097000120826PK
BOX8D: FI:= FI + 1.0; END; GO TO BOX8A; 00098000120826PK
BOX9: GO TO TH[THI]; 00099000
TH1: FOR I:= 1 STEP 1 UNTIL NPP DO 00100000120826PK
X[K,I]:= XP[I]; GO TO BOX10; 00101000120826PK
TH2: FOR I:= 1 STEP 1 UNTIL N DO 00102000120826PK
DEV[R1,I]:= XP[I]; 00103000120826PK
FOR I:= 1 STEP 1 UNTIL P DO 00104000120826PK
COF[R1,I]:= -XP[N+I]; THI:= 3; GO TO TH1; 00105000120826PK
TH3: GO TO BOX11; 00106000
BOX10: GO TO AL[ALI]; 00107000
AL1: OMI:= ALI:= 2; GO TO BOX6; 00108000120826PK
AL2: IF K < M THEN BEGIN K:= K + 1; GO TO BOX4; END 00109000120826PK
ELSE GO TO BOX12; 00110000
BOX11: GO TO NU[NUI]; 00111000
NU1: NUI:= 2; GO TO BOX14; 00112000120826PK
NU2: SS:= DK/DENOM; SSQ:= SS * SS; 00113000120826PK
STD[R1]:= SS; GO TO BOX14; 00114000120826PK
BOX12: GO TO BE[BEI]; 00115000
BE1: FOR I:= 1 STEP 1 UNTIL M DO BEGIN 00116000120826PK
FOR J:= 1 STEP 1 UNTIL P DO 00117000120826PK
A[I,J]:= X[I,N+J] END; 00118000120826PK
GMDT:= 1.0; FOR I:= 1 STEP 1 UNTIL M DO 00119000120826PK
GMDT:= GMDT * (GF[I]/ENF[I]); 00120000120826PK
GMDT:= GMDT * GMDT; DEI:= BEI:= THI:= 2; 00121000120826PK
K:= K + 1; GO TO BOX13; 00122000120826PK
BE2: GO TO BOX11; 00123000
BOX13: GO TO GA[GAI]; 00124000
GA1: GO TO BOX11; 00125000
GA2: FOR I:= 1 STEP 1 UNTIL P DO BEGIN 00126000120826PK
FOR J:= I STEP 1 UNTIL P DO BEGIN 00127000120826PK
SUM:= 0.0; 00128000120826PK
FOR NII:= 1 STEP 1 UNTIL M DO 00129000120826PK
SUM:= SUM + A[NII,I] * A[NII,J]; 00130000120826PK
CV[I,J]:= SUM END END; 00131000120826PK
FOR I:= 1 STEP 1 UNTIL P DO 00132000120826PK
CV[P2,I]:= SQRT(CV[I,I]); GAI:= 1; GO TO BOX11; 00133000120826PK
BOX14: GO TO RH[RHI]; 00134000
RH1: IF RBAR = 0 THEN GO TO FINAL ELSE RBAR:= RBAR - 1; 00135000120826PK
R1:= R1 + 1; THI:= RHI:= 2; GO TO ZE[ZEI]; 00136000120826PK
ZE1: FOR I:= 1 STEP 1 UNTIL N DO 00137000120826PK
X[M2,I]:= Y[R1,I]; 00138000120826PK
FOR I:= 1 STEP 1 UNTIL P DO 00139000120826PK
X[M2,N+I]:= 0.0; GO TO BOX5; 00140000120826PK
ZE2: FOR I:= 1 STEP 1 UNTIL NPP DO 00141000120826PK
X[M2,I]:= Y[R1,I]; GO TO BOX5; 00142000120826PK
RH2: GO TO SI[SII]; 00143000
SI1: GO TO RH1; 00144000
SI2: FOR I:= 1 STEP 1 UNTIL P DO BEGIN 00145000120826PK
FOR J:= I STEP 1 UNTIL P DO 00146000120826PK
VCV[R1,I,J]:= SSQ * CV[I,J] END; 00147000120826PK
FOR I:= 1 STEP 1 UNTIL P DO 00148000120826PK
VCV[R1, P2, I]:= SS * CV[P2,I]; GO TO RH1; 00149000120826PK
FINAL: END ORTHO ; 00150000
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%00151000120429PK
00152000120429PK
PROCEDURE DUMP1 (CAPTION, A); 00800000120826PK
VALUE CAPTION; 00800100120826PK
ALPHA CAPTION; 00800200120826PK
ARRAY A[0]; 00800300120826PK
BEGIN 00800400120826PK
INTEGER 00800500120826PK
UB1, 00800600120826PK
X; 00800700120826PK
00800800120826PK
UB1:= SIZE(A)-1; 00800900120826PK
WRITE (FLUSH, <A6," = ">, CAPTION); 00801000120826PK
WRITE (FLUSH, <6E19.11>, 00801100120826PK
FOR X:= 0 STEP 1 UNTIL UB1 DO A[X]); 00801200120826PK
WRITE (FLUSH); 00801300120826PK
END DUMP1; 00801400120826PK
00810000120826PK
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%00810100120826PK
PROCEDURE DUMP2 (CAPTION, A); 00810200120826PK
VALUE CAPTION; 00810300120826PK
ALPHA CAPTION; 00810400120826PK
ARRAY A[0,0]; 00810500120826PK
BEGIN 00810600120826PK
INTEGER 00810700120826PK
UB1, 00810800120826PK
UB2, 00810900120826PK
X, 00811000120826PK
Y; 00811100120826PK
00811200120826PK
UB1:= SIZE(A)-1; 00811300120826PK
WRITE (FLUSH, <A6," = ">, CAPTION); 00811400120826PK
FOR X:= 0 STEP 1 UNTIL UB1 DO 00811500120826PK
BEGIN 00811600120826PK
UB2:= SIZE(A[X,*])-1; 00811700120826PK
WRITE (FLUSH, <6E19.11>, 00811800120826PK
FOR Y:= 0 STEP 1 UNTIL UB2 DO A[X,Y]); 00811900120826PK
END FOR X; 00812000120826PK
WRITE (FLUSH); 00812100120826PK
END DUMP2; 00812200120826PK
00900000120826PK
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%00900100120826PK
STARTT:= TIME(1); 00901000120826PK
%%READ (PARAMS, /, N, M, P); 00902000120826PK
%%CLOSE(PARAMS); 00903000120826PK
N:= 7; M:= 5; P:= 5; % PARAMETER VALUES ORIGINALLY READ IN 00903100120826PK
R:= 1; AI:= 2; AUI:= 1; 00904000120826PK
ZEI:= 1; MUI:= 2; FN:= N; 00905000120826PK
00906000120429PK
BEGIN %INNER BLOCK 00907000120429PK
LABEL DUMPIT; 00908000120429PK
REAL ARRAY W[1:N], Y[1:R, 1:N+P], Z[1:M, 1:N+P], 00909000120429PK
X[1:M+1, 1:N+P], DEV[1:R, 1:N], COF[1:R, 1:P], STD[1:R], 00910000120429PK
CV[1:P+1, 1:P], VCV[1:R, 1:P+1, 1:P], Q[1:R, 1:M+1], 00911000120429PK
Q2, E, EP [1:R, 1:M], A[1:M, 1:P], GF[1:M+R], ENF[1:M]; 00912000120429PK
%%DUMP FLUSH (N, M, P, R, AI, AUI, ZEI, MUI, 00913000120826PK
%% GMDT, FN, 00914000120826PK
%% W, Y, Z, X, DEV, COF, STD, CV, VCV, Q, Q2, A, GF, ENF, E, 00915000120826PK
%% EP) DUMPIT:1; 00916000120826PK
FOR I:= 1 STEP 1 UNTIL N DO 00917000120826PK
BEGIN 00918000120429PK
W[I]:= 1.0; SUM:= 0.0; 00919000120826PK
FOR J:= 1 STEP 1 UNTIL M DO 00920000120826PK
SUM:= SUM + J*(Z[J,I]:= I*J); 00921000120826PK
Y[1,I]:= SUM; 00922000120826PK
END; 00922500120430PK
ORTHO (W, Y, Z, N, FN, M, P, R, AI, AUI, MUI, ZEI, X, DEV, 00923000120429PK
COF, STD, CV, VCV, GMDT, Q, Q2, E, EP, A, GF, ENF); 00924000120429PK
00924500120430PK
WRITE(FLUSH, <"N = ",I13,/, 00924510120826PK
"M = ",I13,/, 00924520120826PK
"P = ",I13,/, 00924530120826PK
"R = ",I13,/, 00924540120826PK
"AI = ",I13,/, 00924550120826PK
"AUI = ",I13,/, 00924560120826PK
"ZEI = ",I13,/, 00924570120826PK
"MUI = ",I13,/, 00924580120826PK
"GMDT = ",E19.11,/, 00924590120826PK
"FN = ",E19.11,/>, 00924600120826PK
N, M, P, R, AI, AUI, ZEI, MUI, GMDT, FN); 00924610120826PK
DUMP1("W ", W); 00924620120826PK
DUMP2("Y ", Y); 00924630120826PK
DUMP2("Z ", Z); 00924640120826PK
DUMP2("X ", X); 00924650120826PK
DUMP2("DEV ", DEV); 00924660120826PK
DUMP2("COF ", COF); 00924670120826PK
DUMP1("STD ", STD); 00924680120826PK
DUMP2("CV ", CV); 00924690120826PK
DUMP2("VCV ", VCV[1,*,*]); 00924700120826PK
DUMP2("Q ", Q); 00924710120826PK
DUMP2("Q2 ", Q2); 00924720120826PK
DUMP2("A ", A); 00924730120826PK
DUMP1("GF ", GF); 00924740120826PK
DUMP1("ENF ", ENF); 00924750120826PK
DUMP2("E ", E); 00924760120826PK
DUMP2("EP ", EP); 00924770120826PK
DUMPIT: END INNER BLOCK; 00925000120429PK
WRITE(FLUSH, FT, TIME(0), TIME(2)/60, TIME(3)/60, (TIME(1)-STARTT)/60); 00926000120429PK
END. 00999000

View File

@@ -0,0 +1,338 @@
% 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 File

@@ -0,0 +1,344 @@
%%%%%%%%%%%%%%%%%%%%FILTERED MANUALLY BY P.KIMPEL... 2012-08-25 00000100120825PK
% YUSPARA 00336000
% FREE STANDING YU/SPARROW ROUTINE -- 3/5/70 -- P H KIMPEL 00337000
% TEST OF ALGEBRAIC "G" AND "DG" FUNCTIONS -- 8/8/70 -- PHK 00338000
BEGIN 00339000
FILE YUSPAR (KIND=READER, DEPENDENTSPECS); 00340000120825PK
FILE OUTPUT (KIND=PRINTER, MAXRECSIZE=22, FRAMESIZE=48, FILEUSE=OUT); 00341000120825PK
FORMAT 00342000
F1 (//X10"H/R = "F4.1,X15"MU1/MU2 = "I5), 00343000
F2 (2("A[J]"X16"B[J]"X20)/8(2(2(I3,E14.5,X3)X5)X5,2A6"="E13.5/), 00344000
2(I3,E14.5,X3)X55"Q1/Q2"X7"="E13.5), 00345000
F3 (X19"YU/SPARROW ANALYSIS -- ALGEBRAIC """G""" FUNCTIONS"X8,I2, 00346000
2("/"I2)X11"PAGE:"I4), 00347000120825PK
FE (//"RUN TOTALS:"), 00348000
FT (/X9"PROCESSOR TIME: "R8.2,X15"ELAPSED TIME: "R8.2" SEC"); 00349000
INTEGER 00350000
I, J, % INDICIES. 00351000
PT, % PROCESSOR TIME [1/64-S SEC] 00352000
ET, % ELAPSED TIME [SAME]. 00353000
PAG, % SET NUMBER. 00354000
MO, % MONTH-OF-YEAR OF RUN. 00355000
NPOINTS, % TOTAL # POINTS FOR APPROXIMATION. 00356000
NA, % NUMBER OF A-PHASE COEFFICIENTS. 00357000
NB, % NUMBER OF B-PHASE COEFFICIENTS. 00358000
NAB, % NA + NB. 00359000
N1 , % NUMBER OF POINTS ON ARC C1. 00360000
N2, % NUMBER OF POINTS ON ARC C2. 00361000
N12, % NUMBER OF POINTS ON ARC C12. 00362000
DA, YR, % DAY- AND YEAR-OF-RUN. 00362100120825PK
KUSTAR, % COUNTER FOR "USTAR" 00363000
NUSTAR, % COUNT LIMITS FOR USTAR. 00364000
NSUSTAR; % DITTO. 00365000
REAL 00366000
H, % INTERFACE DISTANCE FROM PIPE CENTER. 00370000
R, % PIPE RADIUS (= 1.0) 00371000
MU1, % VISCOSITY OF THE MOST VISCOUS PHASE. 00372000
MU2, % VISCOSITY OF THE LESS VISCOUS PHASE. 00373000
PI, % USUAL MEANING: 3.14159... 00375000
MU12, % MU1/MU2. 00376000
ISUM, % USED BY "INTEGRALUSTAR". 00377000
YSTRT, % DITTO 00378000
YSTEP, % DITTO 00379000
YEND, % DITTO 00380000
XSTEP, % DITTO 00381000
XEND, % DITTO 00382000
XX, % COORDINATE WORKING VARIABLES. 00383000
YY, % DITTO 00384000
USUM, % WORKING VARIABLE VOR "INTEGRALUSTAR" 00385000
Q1Q2, % Q1/Q2. 00386000
X, Y, % WORKING VARIABLES. 00387000
THETA1, % INTERFACE ANGLE FROM VERTICAL DIAMETER. 00388000
ARCC1, % LENGTH OF ARC C1. 00388100
ARCC2, % LENGTH OF ARC C2. 00388200
ARCC12, % LENGTH OF ARC C12. 00388300
DS1, % LENGTH OF ARC BETWEEN POINTS ON C1. 00388400
DS2, % LENGTH OF ARC BETWEEN POINTS ON C2. 00388500
DS12, % LENGTH OF ARC BETWEEN POINTS ON C12. 00388600
TOTALARC; % ARCC1 + ARCC2 + ACRR12. 00388700
LABEL JMP; 00389000
INTEGER ARRAY NDX1, NDX2, NAME1, NAME2 [0:7]; 00391000
ARRAY CALC[0:7]; 00392000
DEFINE Q1Q1FULL = CALC[0]#, Q2Q2FULL = CALC[1]#, Q1STARQ1FULL=CALC[2]#, 00393000
Q2STARQ2FULL = CALC[3]#, TAB22 = CALC[4]#, QMU12 = CALC[5]#, 00394000
Q1Q1STAR = CALC[6]#, Q2Q2STAR = CALC[7]#; 00395000
% BEGIN EXECUTION IN OUTER BLOCK. 00396000
FILL NDX1[*] WITH 1, 3, 4, 7, 8, 11, 12, 15; 00397000
FILL NDX2[*] WITH 2, 5, 6, 9, 10, 13, 14, 17; 00398000
FILL NAME1[*] WITH "Q1/Q1F", "Q2/Q2F", "Q1*/Q1", "Q2*/Q2", "TABLE ", 00399000
"QMU12 ", "Q1/Q1*", "Q2/Q2*"; 00400000
FILL NAME2[*] WITH "ULL ", "ULL ", "FULL ", "FULL ", "2.2 ", 00401000
" ", " ", " "; 00402000
PI := 3.141592653589793; 00403000120825PK
NPOINTS := 100; NA := NB := 17; NAB := NA + NB; 00404000120825PK
R := 1; PAG := 2; 00405000120825PK
DA := TIME(0); YR := DA.[29:6]*10 + DA.[23:6]; 00406000120825PK
DA := DA.[17:6]*100 + DA.[11:6]*10 + DA.[5:6]; 00407000120825PK
FOR I := 31, IF YR MOD 4 = 0 THEN 29 ELSE 28, 31, 00408000120825PK
30, 31, 30, 31, 31, 30, 31, 30, 31 DO 00409000
BEGIN 00410000
MO := MO + 1; 00411000120825PK
IF DA > I THEN DA := DA-I ELSE GO TO JMP; 00412000120825PK
END; 00413000
JMP: 00414000120825PK
FOR H := -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8 DO 00416000120825PK
FOR MU12 := 1, 2, 4, 10, 25, 100, 250, 1000 DO 00417000120825PK
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
COF[1:1, 1:NAB]; 00424000
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; % 100431000120825PK
G := X; % 200432000120825PK
G := Y; % 300433000120825PK
G := X**2 - Y**2; % 400434000120825PK
G := 2.0 * X * Y; % 500435000120825PK
G := X**3 - 3.0 * X * Y**2; % 600436000120825PK
G := 3.0 * X**2 * Y - Y**3; % 700437000120825PK
G := X**4 + Y**4 - 6.0 * X**2 * Y**2; % 800438000120825PK
G := 4.0 * X**3 * Y - 4.0 * X * Y**3; % 900439000120825PK
G := X**5 - 10.0 * X**3 * Y**2 + 5.0 * X * Y**4; %1000440000120825PK
G := Y**5 - 10.0 * Y**3 * X**2 + 5.0 * Y * X**4; %1100441000120825PK
G := X**6 - 15.0 * X**4 * Y**2 + 15.0 * X**2 * Y**4 - Y**6; %1200442000120825PK
G := 6.0 * X**5 * Y + 6.0 * X * Y**5 - 20.0 * X**3 * Y**3; %1300443000120825PK
G := X**7 - 21.0 *X**5 * Y**2 + 35.0 * X**3 * Y**4 - 7.0 * X *%1400444000120825PK
Y**6; 00444100120825PK
G := 7.0 *X**6 * Y - 35.0 * X**4 * Y**3 + 21.0 * X**2 * Y**5 - %1500445000120825PK
Y**7; 00445100120825PK
G := X**8 + Y**8 - 28.0 *X**6 * Y**2 + 70.0 * X**4 * Y**4 - 28.0 00446000120825PK
* X**2 * Y**6; %1600447000120825PK
G := 8.0 * X**7 *Y - 56.0 * X**5 * Y**3 + 56.0 * X**3 * Y**5 00448000120825PK
- 8.0 *X * Y**7; %1700449000120825PK
END G; 00450000
%***********************************************************************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; % 100457000120825PK
DG := 0.0; % 200458000120825PK
DG := 1.0; % 300459000120825PK
DG := -2.0 * Y; % 400460000120825PK
DG := 2.0 * X; % 500461000120825PK
DG := -6.0 * X *Y; % 600462000120825PK
DG := 3.0 * X**2 - 3.0 * Y**2; % 700463000120825PK
DG := 4.0 * Y**3 - 12.0 *X**2 * Y; % 800464000120825PK
DG := 4.0 * X**3 - 12.0 * X * Y**2; % 900465000120825PK
DG := -20.0 *X**3 * Y + 20.0 * X * Y**3; %1000466000120825PK
DG := 5.0 * Y**4 - 30.0 * Y**2 *X**2 + 5.0 * X**4; %1100467000120825PK
DG := -30.0 * X**4 *Y + 60.0 * X**2 * Y**3 - 6.0 * Y**5; %1200468000120825PK
DG := 6.0 * X**5 + 30.0 * X * Y**4 - 60.0 * X**3 * Y**2; %1300469000120825PK
DG := -42.0 * X**5 * Y + 140.0 * X**3 * Y**3 - 42.0 * X * Y**5%1400470000120825PK
; 00470100120825PK
DG := 7.0 * X**6 - 105.0 * X**4 * Y**2 + 105.0 * X**2 * Y**4 00471000120825PK
- 7.0 * Y**6; %1500472000120825PK
DG := 8.0 *Y**7 - 56.0 * X**6 * Y + 280.0 * X**4 * Y**3 00473000120825PK
- 168.0 * X**2 * Y**5; %1600474000120825PK
DG := 8.0 * X**7 - 168.0 * X**5 * Y**2 + 280.0 * X**3 * Y**4 00475000120825PK
- 56.0 * X * Y**6; %1700476000120825PK
END DG; 00477000
%***********************************************************************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; 00483000120825PK
NSUSTAR := IF I = 1 THEN 1 ELSE NA + 1; 00484000120825PK
USUM := 0.0; 00485000120825PK
FOR KUSTAR := NSUSTAR STEP 1 UNTIL NUSTAR DO 00486000120825PK
USUM := USUM + COF[1,KUSTAR] * G(KUSTAR-NSUSTAR+1, X, Y); 00487000120825PK
USTAR := USUM; 00488000120825PK
END USTAR; 00489000
%***********************************************************************00491000
REAL PROCEDURE INTEGRALUSTAR (I); 00492000
VALUE I; INTEGER I; 00493000
BEGIN 00494000
ISUM := 0.0; 00495000120825PK
IF I = 1 THEN BEGIN YSTRT := H; YEND := 0.9999; END 00496000120825PK
ELSE BEGIN YSTRT := -0.9999; YEND := H; END; 00497000120825PK
YSTEP := (YEND - YSTRT)/30.00; 00498000120825PK
FOR YY := YSTRT STEP YSTEP UNTIL YEND DO 00499000120825PK
BEGIN XSTEP := (XX := SQRT(1.0 - YY**2))/15.0; 00500000120825PK
FOR XX := -XX STEP XSTEP UNTIL 0.0 DO 00501000120825PK
ISUM := ISUM + XSTEP * YSTEP * (USTAR (I, XX, YY) 00502000120825PK
- (XX**2 + YY**2)/4); 00503000120825PK
END; 00504000
INTEGRALUSTAR := 2.0 * ISUM; 00505000120825PK
END INTEGRALUSTAR; 00506000
%***********************************************************************00508000
REAL PROCEDURE ARCSIN (X); 00509000
VALUE X; REAL X; 00510000
ARCSIN := ARCTAN(X/SQRT(1.0 - X**2)); 00511000120825PK
%***********************************************************************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 00516000120825PK
+ 1.5 * A *SQRT(1.0 - A**2) + 1.5 * ARCSIN(A)); 00517000120825PK
%***********************************************************************00519000
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
% CODE REDUCED TO ONLY THAT NEEDED TO PERFORM LEAST-SQUARES 00532000
% DETERMINATION OF COEFFICIENTS [4/70]. 00533000
% .... 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; 00543000120825PK
SWITCH TH := TH1, TH2, TH3; SWITCH AL := AL1, AL2; 00544000120825PK
SWITCH OM := OM1, OM2; 00545000120825PK
NPP := N + P; R := 1; M1 := M-1; N2 := N+1; M2 := M+1; 00546000120825PK
R1 := 0; RBAR := R; P2 := P+1; BEI := 1; 00547000120825PK
BOX1: FOR I := 1 STEP 1 UNTIL M DO BEGIN 00548000120825PK
FOR J := 1 STEP 1 UNTIL N DO 00549000120825PK
X[I,J] := Z[I,J] END; 00550000120825PK
BOX2: FOR I := 1 STEP 1 UNTIL M DO BEGIN 00551000120825PK
FOR J := N2 STEP 1 UNTIL NPP DO 00552000120825PK
X[I,J] := 0.0; X[I,N+I] := 1.0 END; GO TO BOX3; 00553000120825PK
BOX3: K := 1; 00554000120825PK
BOX4: THI := 1; 00555000120825PK
BOX5: ALI := OMI := 1; 00556000120825PK
FOR J := 1 STEP 1 UNTIL P DO PK[N+J] := 0.0; 00557000120825PK
BOX6: FOR I := 1 STEP 1 UNTIL N DO 00558000120825PK
PK[I] := X[K,I] * W[I]; 00559000120825PK
BOX7: GO TO OM[OMI]; 00560000
OM1: FOR I := 1 STEP 1 UNTIL K DO BEGIN SUM := 0.0; 00561000120825PK
FOR J := 1 STEP 1 UNTIL NPP DO 00562000120825PK
SUM := SUM + PK[J] * X[I,J]; QK[I] := SUM END; 00563000120825PK
GO TO BOX8; 00564000
OM2: DK2 := 0.0; FOR I := 1 STEP 1 UNTIL NPP DO 00565000120825PK
DK2 := DK2 + PK[I] * X[K,I]; 00566000120825PK
DK := SQRT(DK2); 00567000120825PK
FOR I := 1 STEP 1 UNTIL NPP DO 00568000120825PK
X[K,I] := X[K,I]/DK;; 00569000120825PK
OMI := 1; GO TO BOX6; 00570000120825PK
BOX8: FOR I := 1 STEP 1 UNTIL K-1 DO 00571000120825PK
QK[I] := -QK[I]; QK[K] := 1.0; 00572000120825PK
FOR I := 1 STEP 1 UNTIL NPP DO BEGIN 00573000120825PK
SUM := 0.0; FOR J := 1 STEP 1 UNTIL K DO 00574000120825PK
SUM := SUM + X[J,I] * QK[J]; 00575000120825PK
XP[I] := SUM END; GO TO BOX9; 00576000120825PK
BOX9: GO TO TH[THI]; 00577000
TH1: FOR I := 1 STEP 1 UNTIL NPP DO 00578000120825PK
X[K,I] := XP[I]; GO TO BOX10; 00579000120825PK
TH2: FOR I := 1 STEP 1 UNTIL P DO 00580000120825PK
COF[I] := -XP[N+I]; THI := 3; GO TO TH1; 00581000120825PK
TH3: GO TO BOX14; 00582000
BOX10: GO TO AL[ALI]; 00583000
AL1: OMI := ALI := 2; GO TO BOX6; 00584000120825PK
AL2: IF K < M THEN BEGIN K := K + 1; GO TO BOX4; END 00585000120825PK
ELSE GO TO BOX12; 00586000
BOX12: GO TO BE[BEI]; 00587000
BE1: BEI := THI := 2; K := K + 1; GO TO BOX14; 00588000120825PK
BE2: GO TO BOX14; 00589000
BOX14: IF RBAR = 0 THEN GO TO FINAL ELSE RBAR := RBAR - 1; 00590000120825PK
R1 := R1 + 1; THI := 2; 00591000120825PK
FOR I := 1 STEP 1 UNTIL N DO 00592000120825PK
X[M2,I] := Y[R1,I]; 00593000120825PK
FOR I := 1 STEP 1 UNTIL P DO 00594000120825PK
X[M2,N+I] := 0.0; GO TO BOX5; 00595000120825PK
FINAL: END ORTHO ; 00596000
%***********************************************************************00597000
%***********************************************************************00599000
% BEGIN EXECUTION (INNER BLOCK) %00600000
%***********************************************************************00601000
ET := TIME(1); PT := TIME(2); 00602000120825PK
COMMENT INITALIZE & SET UP DATA FOR -ORTHO-; 00604000
%% CALCULATE ARC LENGTHS OF EACH INTERFACE. 00607000
ARCC12 := 2.0 *SQRT(1.0 - H**2); 00608000120825PK
ARCC1 := 2.0 * (THETA1 := PI/2.0 - ARCTAN(2.0*H/ARCC12)); 00609000120825PK
ARCC2 := 2.0 *PI - ARCC1; 00610000120825PK
%% CALCULATE NO. OF POINTS ON EACH ARC. 00611000
N1 := NPOINTS * ARCC1/(TOTALARC := ARCC1 + ARCC2 + 2*ARCC12); 00612000120825PK
N2 := NPOINTS * ARCC2 / TOTALARC; 00613000120825PK
IF (N12 := NPOINTS - N1 - N2) MOD 2 = 0 THEN 00614000120825PK
N12 := N12 DIV 2 00615000120825PK
ELSE BEGIN 00616000
N12 := N12 DIV 2 + 1; 00617000120825PK
N2 := N2 - 1; 00618000120825PK
END; 00619000
%% CALCULATE THE LENGTH OF LINE SEGMENTS BETWEEN POINTS ON EACH 00620000
%% INTERFACE. 00621000
DS1 := ARCC1/(N1 - 1); 00622000120825PK
DS2 := ARCC2/(N2 - 1); 00623000120825PK
DS12 := ARCC12/(N12 - 1); 00624000120825PK
FOR I := 1 STEP 1 UNTIL NAB DO 00626000120825PK
BEGIN 00627000
FOR J := 1 STEP 1 UNTIL N1 DO %% DATA FOR ARC C1. 00628000120825PK
BEGIN 00629000
W[J] := 1; Y[1,J] := 0.25; 00630000120825PK
XX := SIN(YY := DS1*(J-1)-THETA1); YY := COS(YY); 00631000120825PK
Z[I,J] := IF I LEQ NA THEN G(I, XX, YY) ELSE 0; 00632000120825PK
END; 00633000
FOR J := N1+1 STEP 1 UNTIL N1+N2 DO %% DATA FOR ARC C2. 00634000120825PK
BEGIN 00635000
W[J] := DS2/DS1; Y[1,J] := 0.25; 00636000120825PK
XX := SIN(YY := -THETA1 - DS2*(J-N1-1)); YY := COS(YY); 00637000120825PK
Z[I,J] := IF I LEQ NA THEN 0 ELSE G(I-NA, XX, YY); 00638000120825PK
END; 00639000
FOR J := N1+N2+1 STEP 1 UNTIL NPOINTS-N12 DO %% DATA FOR ARC C12.00640000120825PK
BEGIN 00641000
W[J] := W[J+N12] := DS12/DS1; 00642000120825PK
Y[1,J] := 0.25 * ((XX := DS12*(J-N1-N2-1)-ARCC12/2)**2 + H**2)00643000120825PK
* (1 - MU12); 00644000120825PK
Y[1,J+N12] := 0; 00645000120825PK
IF I LEQ NA THEN 00646000120825PK
BEGIN 00647000
Z[I,J] := G(I, XX, H); 00648000120825PK
Z[I,J+N12] := DG(I, XX, H); 00649000120825PK
END ELSE 00650000
BEGIN 00651000
Z[I,J] := -G(I-NA, XX, H) * MU12; 00652000120825PK
Z[I,J+N12] := -DG(I-NA, XX, H); 00653000120825PK
END; 00654000
END; 00655000
END INITIALIZATION; 00656000
00658000
%%% THE NUMBER SMASHING COMMENCES. 00659000
ORTHO (W, Y, Z, NPOINTS, NAB, NAB, COF[1,*]); %######################00660000
%%% UPON EXIT FROM ORTHO, THE YU/SPARROW COEFFICIENTS, A[J], B[J], 00662000
%%%% ARE CONTAINED IN THE VECTOR COF[1,*]. 00663000
Q1Q1FULL := 8/PI * INTEGRALUSTAR(1); 00665000120825PK
Q2Q2FULL := 8/PI * INTEGRALUSTAR(2); 00666000120825PK
Q1STARQ1FULL := EQUATION23(H); 00667000120825PK
Q2STARQ2FULL := 1 - Q1STARQ1FULL; 00668000120825PK
TAB22 := Q2Q2FULL * PI/8; 00669000120825PK
Q1Q2 := Q1Q1FULL/Q2Q2FULL/MU12; 00670000120825PK
QMU12 := Q1Q2 * MU12; 00670500120825PK
Q1Q1STAR := Q1Q1FULL/Q1STARQ1FULL; Q2Q2STAR := Q2Q2FULL/Q2STARQ2FULL;00671000120825PK
IF (PAG := PAG+1) MOD 3 = 0 THEN 00674000120825PK
BEGIN 00674800
WRITE (OUTPUT[SKIP 1]); 00674900120825PK
WRITE (OUTPUT, F3, MO, DA, YR, PAG DIV 3); 00675000
END; 00675100
WRITE (OUTPUT[SPACE 2], F1, H, MU12); 00676000120825PK
WRITE (OUTPUT, F2, FOR I := 0 STEP 1 UNTIL 7 DO 00677000120825PK
[NDX1[I], COF[1,NDX1[I]], NDX1[I], COF[1,NDX1[I]+NA], NDX2[I], 00678000
COF[1,NDX2[I]], NDX2[I], COF[1,NDX2[I]+NA], NAME1[I], 00679000
NAME2[I], CALC[I]], 16, COF[1,16], 16, COF[1,16], Q1Q2); 00680000
WRITE (OUTPUT, FT, (TIME(2) - PT)/60, (TIME(1) - ET)/60); 00681000
END INNER BLOCK; 00682000
WRITE (OUTPUT[SPACE 2], FE); 00683000120825PK
WRITE (OUTPUT[SPACE 2], *, "TOTALS:"); 00684000120825PK
WRITE (OUTPUT, FT, TIME(2)/60, TIME(1)/60); 00685000
END. 00686000