% 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 IN YUSPAR (1,10); 00340000 FILE OUT OUTPUT 1 "YU/SPAR"/"RESULTS" (2,17); 00341000 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("/"A2)X11"PAGE:"I4), 00347000 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 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 ALPHA DA, YR; % DAY AND YEAR-OF-RUN. 00390000 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; 00403000 NPOINTS ~ 100; NA ~ NB ~ 17; NAB ~ NA + NB; 00404000 R ~ 1; PAG ~ 2; 00405000 DA ~ TIME(0); YR ~ DA.[29:12]; 00406000 DA ~ DA.[17:6]|100 + DA.[11:6]|10 + DA.[5:6]; 00407000 FOR I ~ 31, IF (YR.[11:6]|10 + YR.[5:6])MOD 4 = 0 THEN 29 ELSE 28, 31,00408000 30, 31, 30, 31, 31, 30, 31, 30, 31 DO 00409000 BEGIN 00410000 MO ~ MO + 1; 00411000 IF DA > I THEN DA ~ DA-I ELSE GO TO JMP; 00412000 END; 00413000 JMP: REPLACE POINTER(CALC[0],6)+6 BY DA FOR 2 DIGITS; DA ~ CALC[0]; 00414000 00415000 FOR H ~ -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8 DO 00416000 FOR MU12 ~ 1, 2, 4, 10, 25, 100, 250, 1000 DO 00417000 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; % 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); 00492000 VALUE I; INTEGER I; 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 - (XX*2 + YY*2)/4); 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 % 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; 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 %***********************************************************************00599000 % BEGIN EXECUTION (INNER BLOCK) %00600000 %***********************************************************************00601000 ET ~ TIME(1); PT ~ TIME(2); 00602000 COMMENT INITALIZE & SET UP DATA FOR -ORTHO-; 00604000 %% 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 I ~ 1 STEP 1 UNTIL NAB DO 00626000 BEGIN 00627000 FOR J ~ 1 STEP 1 UNTIL N1 DO %% DATA FOR ARC C1. 00628000 BEGIN 00629000 W[J] ~ 1; Y[1,J] ~ 0.25; 00630000 XX ~ SIN(YY ~ DS1|(J-1)-THETA1); YY ~ COS(YY); 00631000 Z[I,J] ~ IF I { NA THEN G(I, XX, YY) ELSE 0; 00632000 END; 00633000 FOR J ~ N1+1 STEP 1 UNTIL N1+N2 DO %% 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 Z[I,J] ~ IF I { NA THEN 0 ELSE G(I-NA, XX, YY); 00638000 END; 00639000 FOR J ~ N1+N2+1 STEP 1 UNTIL NPOINTS-N12 DO %% 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 - MU12); 00644000 Y[1,J+N12] ~ 0; 00645000 IF I { NA THEN 00646000 BEGIN 00647000 Z[I,J] ~ G(I, XX, H); 00648000 Z[I,J+N12] ~ DG(I, XX, H); 00649000 END ELSE 00650000 BEGIN 00651000 Z[I,J] ~ -G(I-NA, XX, H) | MU12; 00652000 Z[I,J+N12] ~ -DG(I-NA, XX, H); 00653000 END; 00654000 END; 00655000 END INITIALIZATION; 00656000 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/PI | INTEGRALUSTAR(1); 00665000 Q2Q2FULL ~ 8/PI | INTEGRALUSTAR(2); 00666000 Q1STARQ1FULL ~ EQUATION23(H); 00667000 Q2STARQ2FULL ~ 1 - Q1STARQ1FULL; 00668000 TAB22 ~ Q2Q2FULL | PI/8; 00669000 Q1Q2 ~ Q1Q1FULL/Q2Q2FULL/MU12; 00670000 QMU12 ~ Q1Q2 | MU12; 00670500 Q1Q1STAR ~ Q1Q1FULL/Q1STARQ1FULL; Q2Q2STAR ~ Q2Q2FULL/Q2STARQ2FULL;00671000 00672000 00673000 IF (PAG ~ PAG+1) MOD 3 = 0 THEN 00674000 BEGIN 00674800 WRITE (OUTPUT[PAGE]); 00674900 WRITE (OUTPUT, F3, MO, DA, YR, PAG DIV 3); 00675000 END; 00675100 WRITE (OUTPUT[DBL], F1, H, MU12); 00676000 WRITE (OUTPUT, F2, FOR I ~ 0 STEP 1 UNTIL 7 DO 00677000 [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[DBL], FT); 00683000 WRITE (OUTPUT[DBL], *, "TOTALS:"); 00684000 WRITE (OUTPUT, FT, TIME(2)/60, TIME(1)/60); 00685000 END. 00686000