/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:33:13 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv pf=,p_dtgrec s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include #include #include #include #include "p_dtgrec.h" #include "drdtgrec.h" /* File: drDTGrec.[f|for] contains drDTGrec, DTGfcn. */ #include /* PARAMETER translations */ #define MB 16 #define MP 28 #define MT 336 #define MX 5 #define MY 5 #define ZERO 0.0e0 /* end of PARAMETER translations */ int main( ) { char title[76]; LOGICAL32 skip; long int bdry[MB][4], filmod, i, info[3], ip[MP], j, k, kf, kount, ncont, nt, triang[MT]; double delx, dely, dz[MP][2], dzgrd[2][MY][MX], dztrue[MP][2], dzz[2], err, errdz[2], errm, errmdz[2], rms, rmsdz[2], sum1, sum2[2], w[MP], xx, yy, z[MP], zfill, zgrd[MY][MX], zz; static LOGICAL32 prtmax = TRUE; static LOGICAL32 prtt = TRUE; static LOGICAL32 prtr = FALSE; static LOGICAL32 prttpd = FALSE; static char line1[57] = " ------------------------------------------------------"; static char line2[20] = "-------------------"; static double x[MP]={-0.76059e0,-0.02286e0,-0.44790e0,0.15068e0, -0.87287e0,-0.23390e0,0.06093e0,-0.84142e0,-0.69173e0,-0.56613e0, -0.42243e0,-0.17249e0,-0.06484e0,0.48286e0,-0.88784e0,0.10277e0, 0.69087e0,-0.84292e0,0.08784e0,-0.95068e0,0.02496e0,0.94973e0, 0.04588e0,-0.51667e0,-0.77561e0,0.90000e0,-0.70830e0,-0.40500e0}; static double y[MP]={-0.31421e0,0.75657e0,0.14321e0,0.42353e0, -0.62983e0,0.12326e0,0.34054e0,-0.55144e0,0.34158e0,-0.67143e0, 0.43087e0,0.96081e0,-0.68130e0,0.12095e0,-0.10663e0,0.29219e0, 0.31028e0,0.12934e0,0.10709e0,-0.42307e0,0.49895e0,0.68597e0, -0.78215e0,-0.12362e0,-0.88827e0,-0.60000e0,-0.88620e0,0.08600e0}; static double xylim[4]={-1.0e0,1.0e0,-1.0e0,1.0e0}; static long np = MP; static long nx = MX; static long ny = MY; static LOGICAL32 wantpd = TRUE; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Dzz = &dzz[0] - 1; double *const Errdz = &errdz[0] - 1; double *const Errmdz = &errmdz[0] - 1; long *const Info = &info[0] - 1; long *const Ip = &ip[0] - 1; double *const Rmsdz = &rmsdz[0] - 1; double *const Sum2 = &sum2[0] - 1; long *const Triang = &triang[0] - 1; double *const W = &w[0] - 1; double *const X = &x[0] - 1; double *const Xylim = &xylim[0] - 1; double *const Y = &y[0] - 1; double *const Z = &z[0] - 1; /* end of OFFSET VECTORS */ /*>> 2001-07-16 DRDTGREC Krogh Added exponent 0 to some constants. *>> 2001-05-24 DRDTGREC Krogh Minor change for making .f90 version. *>> 2001-05-22 DRDTGREC Krogh Minor change for making .f90 version. *>> 1997-07-01 DRDTGREC Krogh Reversed subscripts in B (CLL suggestion) *>> 1997-06-19 DRDTGREC Krogh Minor changes to get proper C conversion. *>> 1997-06-18 DRDTGREC CLL Added Bdry and MB to arg list of ?TGREC. * Also ZFILL has extended interpretation. *>> 1996-03-04 DRDTGREC CLL *>> 1996-02-02 DRDTGREC CLL *>> 1995-09-26 DRDTGREC CLL Editing for inclusion into MATH77. *>> 1991-11-20 DRDTGREC CLL Minor editing. Adding use of ?SORTP. * C. L. LAWSON, JPL, 1977 MAR 2, CHANGED APR 7 * C.L.L., 1979 MAR 3. MINOR CHANGES. * * PORTABLE DEMONSTRATION DRIVER FOR SUBROUTINES THAT * CONSTRUCT A TRIANGULAR GRID AND THEN DO LOOK-UP AND INTERPOLATION * USING THE GRID. THE INTERPOLATED SURFACE HAS C1 CONTINUITY. * ------------------------------------------------------------------ * The amount of output by this program is controlled by a set of * logical variables. * * PRTMAX Print summary of max and RMS errors over the rectangular * grid. * PRTT Call DTGPRG to print the initial X() and Y() data, the * pointers in TRIANG() defining the triangular grid, and the * pointers in Bdry() identifying the boundary of the triangular * grid. * PRTR Print true and computed values and errors at all nodes of the * rectangular grid. * PRTTPD Print the true and computed values of partial derivatives, * and the errors at all of the initial (x,y) points. * ------------------------------------------------------------------ * THE CALL RELATIONSHIPS ARE AS FOLLOWS.. * * DRDTGREC DEMONSTRATION DRIVER. * DTGGRD CONSTRUCT TRIANGULAR GRID * ?SORTP Sorting. From the JPL Math77 library. * DTGANG COMPUTE PSEUDOANGLES * ?TRSET ** THESE FOUR SUBRS ARE USED TO STORE AND * DTGGET ** FETCH POINTERS IN THE ARRAY TRIANG() * ?TRPUT ** THAT DEFINES THE TRIANGULAR GRID. * ?TRSIZ ** * ?TRADJ TEST AND MODIFY TRIANGULAR GRID * ?TRSET * DTGGET * ?TRPUT * DTGPRG print DESCRIPTION OF TRIANGULAR GRID * DTGGET * DTGFCN PROVIDE VALUES OF TEST functionS. * DTGPD ESTIMATE PARTIAL DERIVATIVES AT GRID POINTS * DTGGET * ?TRMOR FIND MORE PTS FOR PARTIAL DERIV ESTIMATION. * DTGGET * ?TRLS WEIGHT,SCALE,STABILIZE,SOLVE FOR LOCAL L.S. FIT * ?ROTG SETUP FOR GIVENS ROTATION. FROM THE BLAS. * ?ROT APPLY GIVENS ROTATION. FROM THE BLAS. * DTGREC BUILD RECTANGULAR GRID OF INTERPOLATED VALUES. * ?TRFI Lookup and interpolation in triangular grid. * ?TRFND LOOK-UP IN TRIANGULAR GRID * DTGGET * ?TRC0 C0 INTERPOLATION IN A TRIANGLE * ?TRC1 C1 INTERPOLATION IN A TRIANGLE * ------------------------------------------------------------------ * * REMARKS ON DIMENSIONS.. * * IN TERMS OF THE PARAMETERS MP, MT, MB, MX, AND MY, * THE REQUIRED DIMENSIONS ARE.. * double precision ZGRD(MX,MY), DZGRD(MX,MY,2) * double precision X(MP),Y(MP),Z(MP),DZ(2,MP),W(MP),XYLIM(4) * double precision DZZ(2),ERRDZ(2) * double precision DZTRUE(2,MP) * integer IP(MP), TRIANG(MT), Bdry(4,MB) * * (NOTE THAT SOME OF THESE DECLARED VARIABLES ARE ONLY * USED IN THIS DEMO DRIVER ITSELF AND ARE NOT ESSENTIAL * FOR THE GENERAL USE OF THE SUBROUTINES BEING DEMONSTRATED.) * * MP DENOTES THE MAX NO. OF DATA POINTS TO BE HANDLED. * THE MAX NO. OF TRIANGLES CANNOT EXCEED 2*MP. SINCE WE * STORE 6 POINTERS PER TRIANGLE IT FOLLOWS THAT MT CAN BE * SET AS MT = (12*MP) * MB MUST BE LARGER BY 1 THAN THE MAX NO. OF BOUNDARY PTS * IN ANY SUBSET OF THE GIVEN SET OF (X,Y) DATA. THE EXPECTED * NO. OF BNDRY PTS IN A RANDOM SAMPLE OF PTS FROM A UNIFORM * DISTRIBUTION ON DISC IS ABOUT 3.3 * (CUBE ROOT OF MP). * TWICE THIS EXPECTED VALUE PROVIDES AN ADEQUATE BOUND FOR * MOST CASES. SETTING MB = MP+1 WOULD BE ABSOLUTELY SAFE * BUT LARGER THAN NECESSARY IN MOST CASES. * SETTING MB OR MT TOO SMALL WILL BE DETECTED AS ERROR * CONDITIONS BY THE SUBROUTINES. * MX AND MY GIVE THE DIMENSIONS OF THE RECTANGULAR GRID ONTO * WHICH VALUES WILL BE INTERPOLATED UNDER CONTROL OF DTGREC. * ------------------------------------------------------------------ *--D replaces "?": DR?TGREC, ?TGGRD, ?TGPRG, ?TGFCN, ?TGPD, ?TGREC *--& ?TGGET, ?TGANG * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ printf(" Program DRDTGREC. Demo driver for DTGREC,\n DTGGRD, DTGPD, etc.\n\n"); dtggrd( x, y, np, ip, w, triang, MT, bdry, MB, &nt, info ); if (Info[1] != 0) { printf(" Error return from DTGGRD. INFO(1) =%5ld\n", Info[1]); exit(0); } if (prtt) dtgprg( x, y, np, triang, bdry, Info[2], nt ); for (kf = 0; kf <= 4; kf++) { for (i = 1; i <= np; i++) { dtgfcn( kf, X[i], Y[i], title,76, &Z[i], &dztrue[i - 1][0], &dztrue[i - 1][1] ); } printf("\n%s%s\n%s\n%s%s\n", line1, line2, title , line1, line2); dtgpd( x, y, z, dz, np, triang, nt, ip ); if (prttpd) { printf("\n PARTIAL DERIVS ESTIMATED BY DTGPD AT THE GIVEN DATA POINTS\n\n" " I X Y Z DZ1 DZTRUE1 ERR1\n " "DZ2 DZTRUE2 ERR2\n\n"); for (i = 1; i <= np; i++) { for (j = 1; j <= 2; j++) { Errdz[j] = dz[i - 1][j - 1] - dztrue[i - 1][j - 1]; } printf(" %4ld%8.4f%8.4f%8.4f", i, X[i], Y[i], Z[i]); for (j = 1; j <= 2; j++) { printf("%15.6f%15.6f%10.2e", dz[i - 1][j - 1], dztrue[i - 1][j - 1], Errdz[j]); } printf("\n \n"); } } for (ncont = 0; ncont <= 1; ncont++) { for (filmod = 0; filmod <= 1; filmod++) { if (filmod == 0) { zfill = ZERO; } else { zfill = 999.0e0; } dtgrec( x, y, z, dz, np, triang, nt, bdry, MB, xylim, nx, ny, zfill, (double*)zgrd, MX, MY, ncont, wantpd, (double*)dzgrd ); if (prtr || prtmax) { printf("\n\n Using a C%1ld interpolation method.\n", ncont); if (zfill == ZERO) { printf(" Assigning EXTRAPOLATED VALUES to points outside\n " " the convex hull of the given data.\n"); } else { printf(" Assigning the value %7.1f to points outside\n " " the convex hull of the given data.\n", zfill); } if (prtr) { printf("\n VALUES AND PARTIAL DERIVATIVES INTERPOLATED AT\n " " LATTICE POINTS OF A RECTANGULAR GRID.\n\n I J X Y" " Z_INTERP Z_TRUE Z_ERR\n " " DZ1_INTERP DZ1_TRUE DZ1_ERR\n " " DZ2_INTERP DZ2_TRUE DZ2_ERR\n\n"); } delx = (Xylim[2] - Xylim[1])/(nx - 1); dely = (Xylim[4] - Xylim[3])/(ny - 1); xx = Xylim[1]; sum1 = 0.0e0; Sum2[1] = 0.0e0; Sum2[2] = 0.0e0; kount = 0; errm = 0.0e0; Errmdz[1] = 0.0e0; Errmdz[2] = 0.0e0; skip = TRUE; for (i = 1; i <= nx; i++) { yy = Xylim[3]; for (j = 1; j <= ny; j++) { dtgfcn( kf, xx, yy, title,76, &zz, &Dzz[1], &Dzz[2] ); if (zgrd[j - 1][i - 1] != zfill) { err = zgrd[j - 1][i - 1] - zz; if (errm < fabs( err )) errm = fabs( err ); sum1 += SQ(err); kount += 1; for (k = 1; k <= 2; k++) { Errdz[k] = dzgrd[k - 1][j - 1][i - 1] - Dzz[k]; if (Errmdz[k] < fabs( Errdz[k] )) Errmdz[k] = fabs( Errdz[k] ); Sum2[k] += SQ(Errdz[k]); } if (prtr) { printf(" %4ld%4ld%6.2f%6.2f%11.6f%11.6f%10.2e\n %11.6f%11.6f%10.2e\n %11.6f%11.6f%10.2e\n", i, j, xx, yy, zgrd[j - 1][i - 1], zz, err, dzgrd[0][j - 1][i - 1], Dzz[1], Errdz[1], dzgrd[1][j - 1][i - 1], Dzz[2], Errdz[2]); } skip = TRUE; } else { if (prtr) { if (skip) { printf(" \n"); skip = FALSE; } } } yy += dely; } xx += delx; } rms = sqrt( sum1/kount ); for (k = 1; k <= 2; k++) { Rmsdz[k] = sqrt( Sum2[k]/kount ); } if (prtmax) { printf("\n MAX ERR Z MAX ERR ZX MAX ERR ZY\n %17.6g%17.6g%17.6g\n\n " "RMS ERR Z RMS ERR ZX RMS ERR ZY\n %17.6g%17.6g%17.6g\n\n", errm, Errmdz[1], Errmdz[2], rms, Rmsdz[1], Rmsdz[2]); } } } } } exit(0); } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dtgfcn( long kf, double x, double y, char *title,int title_s, double *z, double *zx, double *zy) { static char title1[4-(0)+1][57], title2[4-(0)+1][20]; long int i, _i, _r; static int _aini = 1; #define NCHRTMPS 1 CHRTMP _c[NCHRTMPS]; ini_chrtmp(_c,NCHRTMPS); if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */ { static char* _itmp0[] = {" CONSTANT FUNCTION Z = 2", " "," LINEAR FUNCTION Z = ( 1 + 2*X + 3*Y ) / 6"," ", " QUADRATIC FUNCTION Z = ( -1 + 2*X - 3*Y + 4*X**2 - ", "X*Y + 9*Y**2 ) / 10"," CUBIC FUNCTION Z = ( 9*X**3 - 2*(X**2)*Y + 3*X*Y**2", " - 4 * Y**3 ) / 10"," EXPONENTIAL FUNCTION Z = EXP( -2 * (X**2 + Y**2) )", " "}; for (i = 0, _r = 0; i <= 4; i++) { f_strncpy( title1[i], _itmp0[_r++], 56 ); f_strncpy( title2[i], _itmp0[_r++], 19 ); } } _aini = 0; } /*>> 1995-09-26 CLL Editing for inclusion into MATH77. * C.L.Lawson, JPL, 1976 Dec 10. Edited comments 1979 Mar 3. * This subr evaluates a function and its first partial derivs as * selected by KF. KF can be from 0 to 4. * Input is KF, X, and Y. Output is TITLE, Z, ZX, and ZY. * ------------------------------------------------------------------ *--S replaces "?": ?TRANG * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ f_strncpy( title, f_concat(&_c[0],title1[kf],title2[kf],NULL), title_s - 1 ); if (kf == 0) { /* KF=0 CONSTANT FCN. */ *z = 2.0e0; *zx = 0.0e0; *zy = 0.0e0; } else if (kf == 1) { /* KF=1 LINEAR FCN. */ *z = (1.0e0 + 2.0e0*x + 3.0e0*y)/6.0e0; *zx = 2.0e0/6.0e0; *zy = 3.0e0/6.0e0; } else if (kf == 2) { /* KF=2 QUADRATIC FCN. */ *z = (-1.0e0 + 2.0e0*x - 3.0e0*y + 4.0e0*SQ(x) - x*y + 9.0e0* SQ(y))*0.1e0; *zx = (2.0e0 + 8.0e0*x - y)*0.1e0; *zy = (-3.0e0 - x + 18.0e0*y)*0.1e0; } else if (kf == 3) { /* KF=3 CUBIC FCN. */ *z = (9.0e0*CUBE(x) - 2.0e0*(SQ(x))*y + 3.0e0*x*SQ(y) - 4.0e0* CUBE(y))*0.1e0; *zx = (27.0e0*SQ(x) - 4.0e0*x*y + 3.0e0*SQ(y))*0.1e0; *zy = (-2.0e0*SQ(x) + 6.0e0*x*y - 12.0e0*SQ(y))*0.1e0; } else if (kf == 4) { /* KF=4 EXPONENTIAL FCN. */ *z = exp( -(SQ(x) + SQ(y))*2.0e0 ); /* NOTE THAT THE INFLECTION POINT OF THIS FCN IN ANY * RADIAL DIRECTION FROM THE ORIGIN IS AT R = .5 * */ *zx = -4.0e0*x**z; *zy = -4.0e0*y**z; } rel_chrtmp(_c,NCHRTMPS); return; #undef NCHRTMPS } /* end of function */