diff options
Diffstat (limited to 'sys/gio/ncarutil/conlib/concal.f')
| -rw-r--r-- | sys/gio/ncarutil/conlib/concal.f | 340 |
1 files changed, 340 insertions, 0 deletions
diff --git a/sys/gio/ncarutil/conlib/concal.f b/sys/gio/ncarutil/conlib/concal.f new file mode 100644 index 00000000..e021fa30 --- /dev/null +++ b/sys/gio/ncarutil/conlib/concal.f @@ -0,0 +1,340 @@ + SUBROUTINE CONCAL (XD,YD,ZD,NT,IPT,NL,IPL,PDD,ITI,XII,YII,ZII, + 1 ITPV) +C +C +-----------------------------------------------------------------+ +C | | +C | Copyright (C) 1986 by UCAR | +C | University Corporation for Atmospheric Research | +C | All Rights Reserved | +C | | +C | NCARGRAPHICS Version 1.00 | +C | | +C +-----------------------------------------------------------------+ +C +C +C +C THIS SUBROUTINE PERFORMS PUNCTUAL INTERPOLATION OR EXTRAPO- +C LATION, I.E., DETERMINES THE Z VALUE AT A POINT. +C THE INPUT PARAMETERS ARE +C +C XD,YD,ZD = ARRAYS CONTAINING THE X, Y, AND Z +C COORDINATES OF DATA POINTS, +C NT = NUMBER OF TRIANGLES, +C IPT = INTEGER ARRAY CONTAINING THE POINT NUMBERS OF +C THE VERTEXES OF THE TRIANGLES, +C NL = NUMBER OF BORDER LINE SEGMENTS, +C IPL = INTEGER ARRAY CONTAINING THE POINT NUMBERS OF +C THE END POINTS OF THE BORDER LINE SEGMENTS AND +C THEIR RESPECTIVE TRIANGLE NUMBERS, +C PDD = ARRAY CONTAINING THE PARTIAL DERIVATIVES AT +C THE DATA POINTS, +C ITI = TRIANGLE NUMBER OF THE TRIANGLE IN WHICH LIES +C THE POINT FOR WHICH INTERPOLATION IS TO BE +C PERFORMED, +C XII,YII = X AND Y COORDINATES OF THE POINT FOR WHICH +C INTERPOLATION IS TO BE PERFORMED. +C THE OUTPUT PARAMETER IS +C +C ZII = INTERPOLATED Z VALUE. +C +C DECLARATION STATEMENTS +C +C + DIMENSION XD(1) ,YD(1) ,ZD(1) ,IPT(1) , + 1 IPL(1) ,PDD(1) + DIMENSION X(3) ,Y(3) ,Z(3) ,PD(15) , + 1 ZU(3) ,ZV(3) ,ZUU(3) ,ZUV(3) , + 2 ZVV(3) + REAL LU ,LV + EQUIVALENCE (P5,P50) +C + SAVE +C +C PRELIMINARY PROCESSING +C + IT0 = ITI + NTL = NT+NL + IF (IT0 .LE. NTL) GO TO 100 + IL1 = IT0/NTL + IL2 = IT0-IL1*NTL + IF (IL1 .EQ. IL2) GO TO 150 + GO TO 200 +C +C CALCULATION OF ZII BY INTERPOLATION. +C CHECKS IF THE NECESSARY COEFFICIENTS HAVE BEEN CALCULATED. +C + 100 IF (IT0 .EQ. ITPV) GO TO 140 +C +C LOADS COORDINATE AND PARTIAL DERIVATIVE VALUES AT THE +C IPI 102 VERTEXES. +C IPI 103 +C + JIPT = 3*(IT0-1) + JPD = 0 + DO 120 I=1,3 + JIPT = JIPT+1 + IDP = IPT(JIPT) + X(I) = XD(IDP) + Y(I) = YD(IDP) + Z(I) = ZD(IDP) + JPDD = 5*(IDP-1) + DO 110 KPD=1,5 + JPD = JPD+1 + JPDD = JPDD+1 + PD(JPD) = PDD(JPDD) + 110 CONTINUE + 120 CONTINUE +C +C DETERMINES THE COEFFICIENTS FOR THE COORDINATE SYSTEM +C TRANSFORMATION FROM THE X-Y SYSTEM TO THE U-V SYSTEM +C AND VICE VERSA. +C + X0 = X(1) + Y0 = Y(1) + A = X(2)-X0 + B = X(3)-X0 + C = Y(2)-Y0 + D = Y(3)-Y0 + AD = A*D + BC = B*C + DLT = AD-BC + AP = D/DLT + BP = -B/DLT + CP = -C/DLT + DP = A/DLT +C +C CONVERTS THE PARTIAL DERIVATIVES AT THE VERTEXES OF THE +C TRIANGLE FOR THE U-V COORDINATE SYSTEM. +C + AA = A*A + ACT2 = 2.0*A*C + CC = C*C + AB = A*B + ADBC = AD+BC + CD = C*D + BB = B*B + BDT2 = 2.0*B*D + DD = D*D + DO 130 I=1,3 + JPD = 5*I + ZU(I) = A*PD(JPD-4)+C*PD(JPD-3) + ZV(I) = B*PD(JPD-4)+D*PD(JPD-3) + ZUU(I) = AA*PD(JPD-2)+ACT2*PD(JPD-1)+CC*PD(JPD) + ZUV(I) = AB*PD(JPD-2)+ADBC*PD(JPD-1)+CD*PD(JPD) + ZVV(I) = BB*PD(JPD-2)+BDT2*PD(JPD-1)+DD*PD(JPD) + 130 CONTINUE +C +C CALCULATES THE COEFFICIENTS OF THE POLYNOMIAL. +C + P00 = Z(1) + P10 = ZU(1) + P01 = ZV(1) + P20 = 0.5*ZUU(1) + P11 = ZUV(1) + P02 = 0.5*ZVV(1) + H1 = Z(2)-P00-P10-P20 + H2 = ZU(2)-P10-ZUU(1) + H3 = ZUU(2)-ZUU(1) + P30 = 10.0*H1-4.0*H2+0.5*H3 + P40 = -15.0*H1+7.0*H2-H3 + P50 = 6.0*H1-3.0*H2+0.5*H3 + H1 = Z(3)-P00-P01-P02 + H2 = ZV(3)-P01-ZVV(1) + H3 = ZVV(3)-ZVV(1) + P03 = 10.0*H1-4.0*H2+0.5*H3 + P04 = -15.0*H1+7.0*H2-H3 + P05 = 6.0*H1-3.0*H2+0.5*H3 + LU = SQRT(AA+CC) + LV = SQRT(BB+DD) + THXU = ATAN2(C,A) + THUV = ATAN2(D,B)-THXU + CSUV = COS(THUV) + P41 = 5.0*LV*CSUV/LU*P50 + P14 = 5.0*LU*CSUV/LV*P05 + H1 = ZV(2)-P01-P11-P41 + H2 = ZUV(2)-P11-4.0*P41 + P21 = 3.0*H1-H2 + P31 = -2.0*H1+H2 + H1 = ZU(3)-P10-P11-P14 + H2 = ZUV(3)-P11-4.0*P14 + P12 = 3.0*H1-H2 + P13 = -2.0*H1+H2 + THUS = ATAN2(D-C,B-A)-THXU + THSV = THUV-THUS + AA = SIN(THSV)/LU + BB = -COS(THSV)/LU + CC = SIN(THUS)/LV + DD = COS(THUS)/LV + AC = AA*CC + AD = AA*DD + BC = BB*CC + G1 = AA*AC*(3.0*BC+2.0*AD) + G2 = CC*AC*(3.0*AD+2.0*BC) + H1 = -AA*AA*AA*(5.0*AA*BB*P50+(4.0*BC+AD)*P41)- + 1 CC*CC*CC*(5.0*CC*DD*P05+(4.0*AD+BC)*P14) + H2 = 0.5*ZVV(2)-P02-P12 + H3 = 0.5*ZUU(3)-P20-P21 + P22 = (G1*H2+G2*H3-H1)/(G1+G2) + P32 = H2-P22 + P23 = H3-P22 + ITPV = IT0 +C +C CONVERTS XII AND YII TO U-V SYSTEM. +C + 140 DX = XII-X0 + DY = YII-Y0 + U = AP*DX+BP*DY + V = CP*DX+DP*DY +C +C EVALUATES THE POLYNOMIAL. +C + P0 = P00+V*(P01+V*(P02+V*(P03+V*(P04+V*P05)))) + P1 = P10+V*(P11+V*(P12+V*(P13+V*P14))) + P2 = P20+V*(P21+V*(P22+V*P23)) + P3 = P30+V*(P31+V*P32) + P4 = P40+V*P41 + ZII = P0+U*(P1+U*(P2+U*(P3+U*(P4+U*P5)))) + RETURN +C +C CALCULATION OF ZII BY EXTRATERPOLATION IN THE RECTANGLE. +C CHECKS IF THE NECESSARY COEFFICIENTS HAVE BEEN CALCULATED. +C + 150 IF (IT0 .EQ. ITPV) GO TO 190 +C +C LOADS COORDINATE AND PARTIAL DERIVATIVE VALUES AT THE END +C POINTS OF THE BORDER LINE SEGMENT. +C + JIPL = 3*(IL1-1) + JPD = 0 + DO 170 I=1,2 + JIPL = JIPL+1 + IDP = IPL(JIPL) + X(I) = XD(IDP) + Y(I) = YD(IDP) + Z(I) = ZD(IDP) + JPDD = 5*(IDP-1) + DO 160 KPD=1,5 + JPD = JPD+1 + JPDD = JPDD+1 + PD(JPD) = PDD(JPDD) + 160 CONTINUE + 170 CONTINUE +C +C DETERMINES THE COEFFICIENTS FOR THE COORDINATE SYSTEM +C TRANSFORMATION FROM THE X-Y SYSTEM TO THE U-V SYSTEM +C AND VICE VERSA. +C + X0 = X(1) + Y0 = Y(1) + A = Y(2)-Y(1) + B = X(2)-X(1) + C = -B + D = A + AD = A*D + BC = B*C + DLT = AD-BC + AP = D/DLT + BP = -B/DLT + CP = -BP + DP = AP +C +C CONVERTS THE PARTIAL DERIVATIVES AT THE END POINTS OF THE +C BORDER LINE SEGMENT FOR THE U-V COORDINATE SYSTEM. +C + AA = A*A + ACT2 = 2.0*A*C + CC = C*C + AB = A*B + ADBC = AD+BC + CD = C*D + BB = B*B + BDT2 = 2.0*B*D + DD = D*D + DO 180 I=1,2 + JPD = 5*I + ZU(I) = A*PD(JPD-4)+C*PD(JPD-3) + ZV(I) = B*PD(JPD-4)+D*PD(JPD-3) + ZUU(I) = AA*PD(JPD-2)+ACT2*PD(JPD-1)+CC*PD(JPD) + ZUV(I) = AB*PD(JPD-2)+ADBC*PD(JPD-1)+CD*PD(JPD) + ZVV(I) = BB*PD(JPD-2)+BDT2*PD(JPD-1)+DD*PD(JPD) + 180 CONTINUE +C +C CALCULATES THE COEFFICIENTS OF THE POLYNOMIAL. +C + P00 = Z(1) + P10 = ZU(1) + P01 = ZV(1) + P20 = 0.5*ZUU(1) + P11 = ZUV(1) + P02 = 0.5*ZVV(1) + H1 = Z(2)-P00-P01-P02 + H2 = ZV(2)-P01-ZVV(1) + H3 = ZVV(2)-ZVV(1) + P03 = 10.0*H1-4.0*H2+0.5*H3 + P04 = -15.0*H1+7.0*H2-H3 + P05 = 6.0*H1-3.0*H2+0.5*H3 + H1 = ZU(2)-P10-P11 + H2 = ZUV(2)-P11 + P12 = 3.0*H1-H2 + P13 = -2.0*H1+H2 + P21 = 0.0 + P23 = -ZUU(2)+ZUU(1) + P22 = -1.5*P23 + ITPV = IT0 +C +C CONVERTS XII AND YII TO U-V SYSTEM. +C + 190 DX = XII-X0 + DY = YII-Y0 + U = AP*DX+BP*DY + V = CP*DX+DP*DY +C +C EVALUATES THE POLYNOMIAL. +C + P0 = P00+V*(P01+V*(P02+V*(P03+V*(P04+V*P05)))) + P1 = P10+V*(P11+V*(P12+V*P13)) + P2 = P20+V*(P21+V*(P22+V*P23)) + ZII = P0+U*(P1+U*P2) + RETURN +C +C CALCULATION OF ZII BY EXTRATERPOLATION IN THE TRIANGLE. +C CHECKS IF THE NECESSARY COEFFICIENTS HAVE BEEN CALCULATED. +C + 200 IF (IT0 .EQ. ITPV) GO TO 220 +C +C LOADS COORDINATE AND PARTIAL DERIVATIVE VALUES AT THE VERTEX +C OF THE TRIANGLE. +C + JIPL = 3*IL2-2 + IDP = IPL(JIPL) + X(1) = XD(IDP) + Y(1) = YD(IDP) + Z(1) = ZD(IDP) + JPDD = 5*(IDP-1) + DO 210 KPD=1,5 + JPDD = JPDD+1 + PD(KPD) = PDD(JPDD) + 210 CONTINUE +C +C CALCULATES THE COEFFICIENTS OF THE POLYNOMIAL. +C + P00 = Z(1) + P10 = PD(1) + P01 = PD(2) + P20 = 0.5*PD(3) + P11 = PD(4) + P02 = 0.5*PD(5) + ITPV = IT0 +C +C CONVERTS XII AND YII TO U-V SYSTEM. +C + 220 U = XII-X(1) + V = YII-Y(1) +C +C EVALUATES THE POLYNOMIAL. +C + P0 = P00+V*(P01+V*P02) + P1 = P10+V*P11 + ZII = P0+U*(P1+U*P20) + RETURN + END |