CTITLESAISOV -- ISO-VELOCITY INTERPOLATION OF VELF CRMS FUNCTIONS 00000100 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** C 00000200 CA AUTHOR P. COOPER 00000300 CA DESIGNER P. COOPER 00000400 CA LANGUAGE FORTRAN 00000500 CA SYSTEM IBM AND CARY 00000510 CA WRITTEN 12-14-78 00000600 C REVISED MO-DA-YR BY 00000700 C REVISED 02-24-84 NTS ADD CODE TO RECOMPUTE NUMDO1 AND 00000800 C NUMDO2 ACCORDING TO THE VALUE OF 00000900 C NOSAMP. 00001000 C REVISED 04-19-85 RDK DUAL IBM/CRAY SOURCE CODE. 00001010 C 00001100 CA 00001200 CA 00001300 CA CALL SAISOV (VFT1, VFT2, TIMES1, TIMES2, VELS1, VELS2, 00001400 CA NOSAMP, NUMDO1, NUMDO2) 00001500 CA 00001600 CA INPUT VFT1 = FIRST VELOCITY FUNCTION TRACE R4 00001700 CA INPUT VFT2 = SECOND VELOCITY FUNCTION TRACE R4 00001800 CA IN/OUT TIMES1 = TIMES FOR VELS1 IN VFT1 AND VFT2 I4 00001900 CA IN/OUT TIMES2 = TIMES FOR VELS2 IN VFT1 AND VFT2 I4 00002000 CA INPUT VELS1 = VELOCITIES FROM APPROPRIATE CRMS RECORD R4 00002100 CA INPUT VELS2 = VELOCITIES FROM APPROPRIATE CRMS RECORD R4 00002200 CA INPUT NOSAMP = NUMBER OF SAMPLES IN VFT1, VFT2 I4 00002300 CA IN/OUT NUMDO1 = DIMENSION OF VELS1 I4 00002400 CA IN/OUT NUMDO2 = DIMENSION OF VELS2 I4 00002500 CA 00002600 CA 00002700 CA THIS SUBROUTINE PROJECTS THE INPUT CRMS RECORDS (TIMES1, VELS1,00002800 CA TIMES2, VELS2) FROM ONE VELOCITY FUNCTION INTO A SECOND VELO- 00002900 CA CITY FUNCTION. THIS SUBROUTINE MUST BE RUN BEFORE SAVELI 00003000 CA (WHICH DOES THE ISO-VELOCITY INTERPOLATION) IS RUN. 00003100 CA 00003200 C EJECT 00003300 C 00003400 C THE FIRST COLUMN OF TIMES1 CONTAINS THE TIMES OF THE CRMS 00003500 C RECORD FOR THE FIRST VELOCITY FUNCTION. VELS1 CONTAINS THE 00003600 C VELOCITIES FROM THE SAME RECORD. THE SECOND COLUMN OF TIMES2 00003700 C CONTAINS THE TIMES OF THE CRMS RECORD FOR THE SECOND VELOCITY 00003800 C FUNCTION. VELS2 CONTAINS THE VELOCITIES FROM THE SAME RECORD. 00003900 C 00004000 C THE SECOND COLUMN OF TIMES1 AND THE FIRST COLUMN OF TIMES2 00004100 C ARE COMPUTED OR FOUND IN THIS PROGRAM. THE SECOND COLUMN OF 00004200 C TIMES1 CONTAINS THE TIMES WHERE THE VELOCITIES GIVEN IN VELS1 00004300 C WERE MATCHED IN VFT2 (THE SECOND VELOCITY FUNCTION). THE 00004400 C FIRST COLUMN OF TIMES2 CONTAINS THE TIMES WHERE THE VELOCITIES 00004500 C GIVEN IN VELS2 WERE MATCHED IN VFT1 (THE FIRST VELOCITY 00004600 C FUNCTION). 00004700 C 00004800 C TO MATCH THE CRMS VELOCITIES TO THE VELOCITY FUNCTIONS, A 00004900 C SEARCH OF THE VELOCITY FUNCTION IS MADE. THE SEARCH STARTS 00005000 C FROM THE LAST VELOCITY MATCH (OR THE START OF THE FUNCTION) 00005100 C AND CONTINUES UNTIL THE VELOCITY FUNCTION VELOCITY IS GREATER 00005200 C THAN OR EQUAL TO THE CRMS VELOCITY. THE TIME WHERE THE SEARCH 00005300 C STOPS IS THEN PUT IN THE APPROPRIATE POSITION OF TIMES1 OR 00005400 C TIMES2. THIS POSITION IS DETERMINED BY THE POSITION OF 00005500 C THE VELOCITY IN VELS1 OR VELS2 FOR WHICH WE ARE SEARCHING. 00005600 C 00005700 C THERE ARE THREE SPECIAL CASES. ONE IS WHERE THE VELOCITY 00005800 C BEING SEARCHED FOR IS LESS THAN THE FIRST VELOCITY OF THE 00005900 C VELOCITY FUNCTION. ONE IS WHERE THE VELOCITY BEING SEARCHED 00006000 C FOR IS GREATER THAN THE LAST VELOCITY OF THE VELOCITY FUNCTION.00006100 C AND ONE IS WHERE THERE IS AN INVERSION. THE FIRST TWO CASES 00006200 C ARE HANDLED IN SAVELI. WHEN THERE IS AN INVERSION, THE TIMES 00006300 C WHERE THE INVERSION BEGINS AND ENDS ON EACH VELOCITY FUNCTION 00006400 C ARE RECORDED. MATCHING TIMES FOR ALL POINTS WITHIN THE INVER- 00006500 C SION ARE COMPUTED USING A RELATIVE DISTANCE FROM THE BEGINNING 00006600 C OF THE INVERSION. CHECKS ARE MADE TO BE SURE THAT NO ISO - 00006700 C VELOCITY LINES CROSS THE INVERSION. IF THE INVERSION OCCURS 00006800 C IN ONLY ONE OF THE FUNCTIONS THEN IT IS GRADUALLY TAPERED 00006900 C INTO THE OTHER FUNCTION. AN INVERSION IS MARKED BY A NEGATIVE 00007000 C TIME IN THE FIRST COLUMN. 00007100 C 00007200 C THESE TIMES THAT HAVE BEEN COMPUTED ARE PASSED BACK TO SDNMOC 00007300 C WHERE THEY ARE USED BY SAVELI. 00007400 C 00007500 C EJECT 00007600 SUBROUTINE SAISOV (VFT1, VFT2, TIMES1, TIMES2, VELS1, VELS2, 00007700 * NOSAMP, NUMDO1, NUMDO2) 00007800 C 00007900 IMPLICIT INTEGER (A-Z) 00008000 C 00008100 C REAL ARRAYS IN ARGUMENT LIST 00008200 C 00008300 REAL VFT1 (1) 00008400 REAL VFT2 (1) 00008500 REAL VELS1 (114) 00008600 REAL VELS2 (114) 00008700 C 00008800 C INTEGER ARRAYS IN ARGUMENT LIST AND LOCAL 00008900 C 00009000 INTEGER TIMES1 (114,2) 00009100 INTEGER TIMES2 (114,2) 00009200 INTEGER STRINV (10,2) 00009300 C 00009400 C REAL VARIABLES -- LOCAL 00009500 C 00009600 REAL TEMP3 00009700 REAL SLOPE 00009800 C 00009900 C 00010000 C ADD CODE HERE TO COMPARE NOSAMP WITH MAX-TIME ON TIME-VELOCITY 00010100 C PAIRS. IF MAX-TIME IS LARGER THAN NOSAMP THEN RECOMPUTE 00010200 C NUMDO1 AND NUMDO2 IF NECESSARY 00010300 C 00010400 IF ( NOSAMP+1 .GE. TIMES1(NUMDO1,1) ) GO TO 5 00010500 C 00010600 C FIND THE CORRECT VALUE OF NUMDO1 HERE 00010700 C 00010800 DO 1 JJ=1, NUMDO1 00010900 KK = NUMDO1 + 1 - JJ 00011000 IF ( NOSAMP .LE. TIMES1(KK-1,1 )) GO TO 1 00011100 C 00011200 C CALCIULATE NOW NEW CONTROL POINT RMS-VELOCITY 00011300 C 00011400 SLOPE =(VELS1(KK) - VELS1( KK - 1))/ 00011500 * (TIMES1(KK,1) - TIMES1( KK - 1,1 )) 00011600 VELS1 ( KK) = VELS1 ( KK - 1) 00011700 * + SLOPE*(NOSAMP+1 - TIMES1(KK-1,1)) 00011800 TIMES1(KK,1) = NOSAMP + 1 00011900 SAVEKK = KK 00012000 GO TO 3 00012100 1 CONTINUE 00012200 3 NUMDO1 = SAVEKK 00012300 C 00012400 5 CONTINUE 00012500 IF ( NOSAMP + 1.GE. TIMES2( NUMDO2, 2 )) GO TO 9 00012600 C 00012700 C FIND THE CORRECT VALUE OF NUMDO2 HERE 00012800 C 00012900 DO 6 JJ=1, NUMDO2 00013000 KK = NUMDO2 + 1 - JJ 00013100 IF ( NOSAMP .LE. TIMES2 ( KK - 1,2)) GO TO 6 00013200 C 00013300 C CALCIULATE NOW NEW CONTROL POINT RMS-VELOCITY 00013400 C 00013500 SLOPE =(VELS2 (KK ) - VELS2 ( KK - 1))/ 00013600 * (TIMES2(KK,2) - TIMES2( KK - 1,2 )) 00013700 VELS2 ( KK ) = VELS2 ( KK - 1) 00013800 * + SLOPE*(NOSAMP+1 - TIMES2(KK-1,2)) 00013900 TIMES2(KK,2) = NOSAMP + 1 00014000 SAVEKK = KK 00014100 GO TO 8 00014200 6 CONTINUE 00014300 8 NUMDO2 = SAVEKK 00014400 C 00014500 9 CONTINUE 00014600 C 00014700 C INITIALIZE START TIME 00014800 C DETERMINE WHICH VELOCITY TO START MATCHING 00014900 C 00015000 ENDDO = NOSAMP + 1 00015100 ST = 1 00015200 10 IF (VELS1(ST) .GE. VFT2(1)) GO TO 20 00015300 ST = ST + 1 00015400 GO TO 10 00015500 C 00015600 C INITIALIZE VARIABLES 00015700 C 00015800 20 INDX1 = 1 00015900 INVFLG = 0 00016000 INVIDX = 0 00016100 INVFL1 = 0 00016200 C 00016300 C LOOP FOR MATCHING OF VELS1 VELOCITIES TO VFT2. 00016400 C CHECK FOR INVERSION, IF INVERSION THEN SET 00016500 C TIME OF VELS1 NEGATIVE 00016600 C 00016700 DO 80 00016800 * I = ST, NUMDO1 00016900 IF (I-1 .EQ. 0) GO TO 50 00017000 IF (VELS1(I) .GE. VELS1(I-1)) GO TO 30 00017100 IF (INVFLG .EQ. 0) INVFLG = I-1 00017200 30 IF (INVFLG .EQ. 0) GO TO 50 00017300 IF (VELS1(I) .GE. VELS1(INVFLG)) GO TO 40 00017400 TIMES1(I,1) = -TIMES1(I,1) 00017500 GO TO 80 00017600 C 00017700 C STORE STARTING AND ENDING POINTS OF INVERSIONS 00017800 C 00017900 40 INVIDX = INVIDX + 1 00018000 STRINV(INVIDX,1) = INVFLG 00018100 STRINV(INVIDX,2) = I 00018200 INVFL1 = 1 00018300 INVFLG = 0 00018400 C 00018500 C LOOP TO FIND TIME WHERE VELS1 AND VFT2 MATCH 00018600 C 00018700 50 DO 60 00018800 * J = INDX1, ENDDO 00018900 IF (VELS1(I) .GT. VFT2(J)) GO TO 60 00019000 TIMES1(I,2) = J 00019100 IF (I-1 .EQ. 0) GO TO 70 00019200 IF (VELS1(I) .EQ. VELS1(I-1) .AND. 00019300 * VELS1(I) .EQ. VFT2(J)) GO TO 60 00019400 IF (VELS1(I) .EQ. VELS1(I-1)) TIMES1(I,2) = J - 1 00019500 GO TO 70 00019600 60 CONTINUE 00019700 C 00019800 C IF NO VELOCITY MATCH FOUND GO TO NEXT VELOCITY 00019900 C 00020000 GO TO 80 00020100 C 00020200 C SET START TIME OF NEXT SEARCH 00020300 C 00020400 70 INDX1 = TIMES1(I,2) + 1 00020500 C 00020600 80 CONTINUE 00020700 C 00020800 C LOOP TO MATCH VELS2 VELOCITIES TO VFT1 00020900 C 00021000 ST = 1 00021100 90 IF (VELS2(ST) .GE. VFT1(1)) GO TO 100 00021200 ST = ST + 1 00021300 GO TO 90 00021400 C 00021500 C INITIALIZE VARIABLES 00021600 C 00021700 100 INDX1 = 1 00021800 INVFLG = 0 00021900 IVINDX = 0 00022000 C 00022100 C MATCH VELS2 VELOCITIES TO VFT1, CHECK FOR INVERSION 00022200 C 00022300 DO 310 00022400 * I = ST, NUMDO2 00022500 IF (I-1 .EQ. 0) GO TO 120 00022600 IF (VELS2(I) .GE. VELS2(I-1)) GO TO 110 00022700 IF (INVFLG .EQ. 0) INVFLG = I - 1 00022800 110 IF (INVFLG .EQ. 0) GO TO 120 00022900 IF (VELS2(I) .GE. VELS2(INVFLG)) GO TO 120 00023000 GO TO 310 00023100 C 00023200 C FIND TIME OF VELS2 IN VFT1 00023300 C 00023400 120 DO 130 00023500 * J = INDX1, ENDDO 00023600 IF (VELS2(I) .GT. VFT1(J)) GO TO 130 00023700 TIMES2(I,1) = J 00023800 IF (I-1 .EQ. 0) GO TO 140 00023900 IF (VELS2(I) .EQ. VELS2(I-1) .AND. 00024000 * VELS2(I) .EQ. VFT1(J)) GO TO 130 00024100 IF (VELS2(I) .EQ. VELS2(I-1)) TIMES2(I,1) = J - 1 00024200 GO TO 140 00024300 130 CONTINUE 00024400 C 00024500 C IF NO MATCH FOUND DO NEXT VELOCITY 00024600 C 00024700 GO TO 310 00024800 C 00024900 C CALCULATE TIME TO MATCH IF WAS AN INVERSION 00025000 C 00025100 140 INDX1 = TIMES2(I,1) + 1 00025200 IF (INVFLG .EQ. 0) GO TO 310 00025300 IF (INVFL1 .NE.0) GO TO 200 00025400 TEMP1 = TIMES2(I,1) - TIMES2(INVFLG,1) 00025500 TEMP2 = TIMES2(I,2) - TIMES2(INVFLG,2) 00025600 K = INVFLG 00025700 C 00025800 150 K = K + 1 00025900 IF (K .EQ. I) GO TO 160 00026000 TEMP3 = TIMES2(K,2) - TIMES2(INVFLG,2) 00026100 TEMP3 = TEMP3 / TEMP2 00026200 TIMES2(K,1) = -(TEMP3 * TEMP1 + TIMES2(INVFLG,1)) 00026300 GO TO 150 00026400 C 00026500 C CHECK IF TIMES CROSS INVERSION 00026600 C 00026700 160 TIME1 = IABS(TIMES2(INVFLG+1,1)) 00026800 TIME2 = TIMES2(INVFLG+1,2) 00026900 TIME3 = IABS(TIMES2(I-1,1)) 00027000 TIME4 = TIMES2(I-1,2) 00027100 C 00027200 DO 190 00027300 * L = 1, NUMDO1 00027400 IF (TIMES1(L,2) .EQ. 0) GO TO 190 00027500 IF (TIMES1(L,1) .GT. TIME1 .OR. 00027600 * TIMES1(L,2) .LE. TIME2) GO TO 170 00027700 TIMES1(L,2) = TIME2 00027800 TIMES1(L,1) = -TIMES1(L,1) 00027900 GO TO 190 00028000 C 00028100 170 IF (TIMES1(L,1) .LE. TIME1 .OR. 00028200 * TIMES1(L,1) .GT. TIME3) GO TO 180 00028300 TEMP1 = TIME3 - TIME1 00028400 TEMP2 = TIME4 - TIME2 00028500 TEMP3 = TIMES1(L,1) - TIME1 00028600 TEMP3 = TEMP3 / TEMP1 00028700 TIMES1(L,2) = TEMP3 * TEMP2 + TIME2 00028800 TIMES1(L,1) = -TIMES1(L,1) 00028900 GO TO 190 00029000 C 00029100 180 IF (TIMES1(L,1) .LT. TIME3 .OR. 00029200 * TIMES1(L,2) .GE. TIME4) GO TO 190 00029300 TIMES1(L,2) = TIME4 00029400 TIMES1(L,1) = -TIMES1(L,1) 00029500 C 00029600 190 CONTINUE 00029700 C 00029800 GO TO 300 00029900 C 00030000 C MUST CHECK IF FIRST VFT HAD AN INVERSION AND MAKE 00030100 C SURE NO ISO-VELOCITY LINES CROSS THE INVERSION 00030200 C 00030300 200 IVINDX = IVINDX + 1 00030400 INVER1 = STRINV(IVINDX,1) 00030500 INVER2 = INVFLG 00030600 C 00030700 C CHECK IF FIRST VFT CROSSES INVERSION -- DOWN 00030800 C 00030900 210 IF (TIMES1(INVER1,2) .LE. TIMES2(INVFLG,2)) GO TO 22000031000 TIMES1(INVER1,2) = TIMES2(INVFLG,2) 00031100 TIMES1(INVER1,1) = -TIMES1(INVER1,1) 00031200 INVER1 = INVER1 - 1 00031300 IF (INVER1 .EQ. 0) GO TO 220 00031400 GO TO 210 00031500 C 00031600 C CHECK IF SECOND VFT CROSSES INVERSION -- DOWN 00031700 C 00031800 220 IF (TIMES2(INVER2,1) .LE. 00031900 * IABS(TIMES1(STRINV(IVINDX,1),1))) GO TO 230 00032000 TIMES2(INVER2,1) = -IABS(TIMES1(STRINV(IVINDX,1),1)) 00032100 INVER2 = INVER2 - 1 00032200 IF(INVER2 .EQ. 0) GO TO 230 00032300 GO TO 220 00032400 C 00032500 C SET TOP END TIMES OF INVERSION 00032600 C 00032700 230 TIME1 = IABS(TIMES1(STRINV(IVINDX,1),1)) 00032800 TIME2 = TIMES2(INVFLG,2) 00032900 INVER1 = STRINV(IVINDX,2) 00033000 INVER2 = I 00033100 C 00033200 C CHECK IF FIRST VFT CROSSES INVERSION -- UP 00033300 C 00033400 240 IF (TIMES1(INVER1,2) .GE. TIMES2(I,2)) GO TO 250 00033500 IF (TIMES1(INVER1,2) .EQ. 0) GO TO 250 00033600 TIMES1(INVER1,2) = TIMES2(I,2) 00033700 TIMES1(INVER1,1) = -TIMES1(INVER1,1) 00033800 INVER1 = INVER1 + 1 00033900 IF (INVER1 .EQ. NUMDO1) GO TO 250 00034000 GO TO 240 00034100 C 00034200 C CHECK IF SECOND CROSSES INVERSION -- UP 00034300 C 00034400 250 IF (TIMES2(INVER2,1) .GE. 00034500 * IABS(TIMES1(STRINV(IVINDX,2),1))) GO TO 260 00034600 IF (TIMES2(INVER2,1) .EQ. 0) GO TO 260 00034700 TIMES2(INVER2,1) = -IABS(TIMES1(STRINV(IVINDX,2),1)) 00034800 INVER2 = INVER2 + 1 00034900 IF (INVER2 .EQ. NUMDO2) GO TO 260 00035000 GO TO 250 00035100 C 00035200 C SET END TIMES OF INVERSION AND CALCULATE WIDTH 00035300 C OF INVERSION ON EACH VELOCITY 00035400 C 00035500 260 TIME3 = IABS(TIMES1(STRINV(IVINDX,2),1)) 00035600 TIME4 = TIMES2(I,2) 00035700 TEMP1 = TIME3 - TIME1 00035800 TEMP2 = TIME4 - TIME2 00035900 END1 = STRINV(IVINDX,2) 00036000 K = STRINV(IVINDX,1) 00036100 C 00036200 C CALCULATE PROJECTED TIME OF FIRST VFT INTO SECOND 00036300 C 00036400 270 K = K + 1 00036500 IF (K .EQ. END1) GO TO 280 00036600 TEMP3 = IABS(TIMES1(K,1)) - TIME1 00036700 TEMP3 = TEMP3 / TEMP1 00036800 TIMES1(K,2) = TEMP3 * TEMP2 + TIME2 00036900 GO TO 270 00037000 C 00037100 C CALCULATE PROJECTED TIME OF SECOND VFT IN FIRST 00037200 C 00037300 280 K = INVFLG 00037400 290 K = K + 1 00037500 IF (K .EQ. I) GO TO 300 00037600 TEMP3 = TIMES2(K,2) - TIME2 00037700 TEMP3 = TEMP3 / TEMP2 00037800 TIMES2(K,1) = -(TEMP3 * TEMP1 + TIME1) 00037900 GO TO 290 00038000 C 00038100 C RESET INVERSION FLAGS 00038200 C 00038300 300 INVFLG = 0 00038400 IF (IVINDX .EQ. INVIDX) INVFL1 = 0 00038500 C 00038600 310 CONTINUE 00038700 C 00038800 C CHECK IF STILL AN INVERSION IN FIRST VFT 00038900 C 00039000 320 IF (INVFL1 .EQ. 0) GO TO 380 00039100 IVINDX = IVINDX + 1 00039200 INDX1 = STRINV(IVINDX,1) 00039300 INDX2 = STRINV(IVINDX,2) 00039400 C 00039500 C CALCULATE WIDTH OF INVERSION 00039600 C 00039700 TEMP1 = TIMES1(INDX2,1) - TIMES1(INDX1,1) 00039800 TEMP2 = TIMES1(INDX2,2) - TIMES1(INDX1,2) 00039900 K = INDX1 00040000 C 00040100 C PROJECT TIME OF INVERSION 00040200 C 00040300 330 K = K + 1 00040400 IF (K .EQ. INDX2) GO TO 340 00040500 TEMP3 = IABS(TIMES1(K,1)) - TIMES1(INDX1,1) 00040600 TEMP3 = TEMP3 / TEMP1 00040700 TIMES1(K,2) = TEMP3 * TEMP2 + TIMES1(INDX1,2) 00040800 GO TO 330 00040900 C 00041000 C CHECK IF VELOCITIES CROSS INVERSION 00041100 C 00041200 340 TIME1 = IABS(TIMES1(INDX1+1,1)) 00041300 TIME2 = TIMES1(INDX1+1,2) 00041400 TIME3 = IABS(TIMES1(INDX2-1,1)) 00041500 TIME4 = TIMES1(INDX2-1,2) 00041600 C 00041700 DO 370 00041800 * L =1,NUMDO2 00041900 C 00042000 IF (TIMES2(L,2) .GT. TIME2 .OR. 00042100 * TIMES2(L,1) .LE. TIME1) GO TO 350 00042200 TIMES2(L,1) = -TIME1 00042300 GO TO 370 00042400 C 00042500 350 IF (TIMES2(L,2) .LE. TIME2 .OR. 00042600 * TIMES2(L,2) .GT. TIME4) GO TO 360 00042700 TEMP1 = TIME3 - TIME1 00042800 TEMP2 = TIME4 - TIME2 00042900 TEMP3 = TIMES2(L,2) - TIME2 00043000 TEMP3 = TEMP3 / TEMP2 00043100 TIMES2(L,1) = -(TEMP3 * TEMP1 + TIME1) 00043200 GO TO 370 00043300 C 00043400 360 IF (TIMES2(L,2) .LT. TIME4 .OR. 00043500 * TIMES2(L,1) .GE. TIME3) GO TO 370 00043600 IF (TIMES2(L,1) .EQ. 0) GO TO 370 00043700 TIMES2(L,1) = -TIME3 00043800 C 00043900 370 CONTINUE 00044000 C 00044100 C 00044200 C CHECK IF STILL ANOTHER INVERSION IN FIRST VFT 00044300 C 00044400 IF (IVINDX .NE. INVIDX) GO TO 320 00044500 C 00044600 380 CONTINUE 00044700 C 00044710 RETURN 00044800 C 00044810 END 00044900