/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:54 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dtgpd.h" #include #include /* PARAMETER translations */ #define MFMAX 16 #define MW (MFMAX + 5) /* end of PARAMETER translations */ /* COMMON translations */ struct t_dtgcm1 { double x0, y0, z0, dsq[MFMAX]; long int avail, mfmin, mfit, firstl, lastl, ke[4][MFMAX], jused[MFMAX]; LOGICAL32 newpt; } dtgcm1; /* end of COMMON translations */ void /*FUNCTION*/ dtgpd( double x[], double y[], double z[], double dz[][2], long np, long triang[], long nt, long iwork[]) { LOGICAL32 hitbdy; long int i, ip1, ipcent, it, iv1, iv2, ivert, j1, j2, jj1, jj2, jp, kp, l, limit, lolim, tnow, tri[6], trsave[6], _i, _r; static long int add1[3], sub1[3]; double w[6][21], x1, x2, y1, y2; static int _aini = 1; /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Add1 = &add1[0] - 1; double *const Dsq = &dtgcm1.dsq[0] - 1; long *const Iwork = &iwork[0] - 1; long *const Jused = &dtgcm1.jused[0] - 1; long *const Sub1 = &sub1[0] - 1; long *const Tri = &tri[0] - 1; long *const Trsave = &trsave[0] - 1; double *const X = &x[0] - 1; double *const Y = &y[0] - 1; double *const Z = &z[0] - 1; /* end of OFFSET VECTORS */ if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */ Add1[1] = 2; Add1[2] = 3; Add1[3] = 1; Sub1[1] = 3; Sub1[2] = 1; Sub1[3] = 2; _aini = 0; } /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. * File: DTGPD.[f|for] contains DTGPD, DTGMOR, DTGLS. *>> 2007-02-28 DTGPD Krogh Massive changed to remove assigned go to's. *>> 2005-12-07 DTGPD Krogh Removed unused label. *>> 1996-05-10 DTGPD Krogh MW removed from dim. declarator in C version *>> 1996-03-30 DTGPD Krogh MIN0 => MIN *>> 1996-02-02 DTGPD CLL *>> 1996-01-15 DTGPD CLL *>> 1996-01-11 DTGPD CLL *>> 1995-09-26 DTGPD CLL Editing for inclusion into MATH77. * * THIS SUBR ESTIMATES FIRST PARTIAL DERIVS AT THE GIVEN DATA POINTS. * C.L.LAWSON, JPL, 1976 DEC 21. EDITED COMMENTS 1979 MAR 5. * * --------------------------------------------------------------------- * Method * * THE PARTIAL DERIVS AT POINT JP WILL BE ESTIMATED BY FITTING * A QUADRATIC IN X AND Y TO A SET OF MFIT+1 POINTS CONSISTING OF * POINT JP AND MFIT NEARBY POINTS IN THE TRIANGULAR GRID. * THE VALUE OF MFIT WILL IN GENERAL BE DIFFERENT AT DIFFERENT * POINTS. MFIT WILL BE IN THE RANGE FROM MFMIN TO MFMAX. LET NNEB * DENOTE THE NO. OF IMMEDIATE NEIGHBORING POINTS TO POINT JB. IF * MFMIN .LE. NNEB .LE. MFMAX THEN MFIT=NNEB. * IF NNEB .LT. MFMIN THEN MFIT = MFMIN. * IF NNEB .GT. MFMAX THEN MFIT = MFMAX. * --------------------------------------------------------------------- * Subroutine Arguments * * X(), Y(), Z() [in] * DZ(2,) [out] * NP [in] * TRIANG(), NT [in] * (IWORK(I),I=1,NP) [scratch] INITIALLY SET TO ZERO BY DTGPD. * DTGPD WILL SET IWORK(JP)= 1 WHEN IT ESTIMATES THE PARTIAL * DERIVATIVES AT POINT JP. * * --------------------------------------------------------------------- * COMMON BLOCK * * /DTGCM1/ provides communication between DTGPD and DTGMOR. * --------------------------------------------------------------------- * Internal variables * * AVAIL POINTS TO AVAILABLE SPACE IN KE(). * MFMAX Dimensioning parameter used in the following declarations: * INTEGER KE(MFMAX ,4), JUSED(MFMAX) * DOUBLE PRECISION W(MFMAX+5,6), DSQ(MFMAX) * MFMIN LOWER BOUND FOR MFIT. SET TO A CONSTANT VALUE BY DTGPD. * MFIT NO. OF NEIGHBORING POINTS USED WITH EACH POINT JP * TO ESTIMATE PARTIALS AT POINT JP. GENERALLY DIFFERS * FOR EACH POINT JP. * FIRSTL,LASTL POINTERS TO FIRST AND LAST ELEMENTS OF THE * LIST CONTAINED IN KE(). * KE(,) IS A DOUBLY LINKED LIST DESCRIBING A CYCLE OF EDGES * SURROUNDING POINT JP. AN EDGE IS DEFINED BY ITS TWO ADJACENT * TRIANGLES. * KE(L,1) FWD POINTER (COUNTERCLOCKWISE). * KE(L,2) BKWD POINTER (CLOCKWISE). * KE(L,3) TRIANGLE INSIDE EDGE L. * KE(L,4) TRIANGLE OUTSIDE EDGE L, OR ELSE ZERO IF EDGE * ON THE BOUNDARY. * JUSED() AN UNORDERED SET OF INDICES OF ROINTS USED IN THE LOCAL * FIT THAT ARE NOT IMMEDIATE NEIGHBORS OF POINT JP. * X0,Y0,Z0 COORDINATES OF POINT JP. * DSQ(L) SQUARED DISTANCE FROM POINT JP TO THE MIDPOINT OF SIDE L * NEWPT LOGICAL FLAG. FOR A NEW POINT JP DTGPD SETS * NEWPT TRUE. THIS CAUSES DTGLS TO TRIANGULARIZE THE * LEAST SQUARES SYSTEM FOR JP. DTGLS THEN TESTS THE * CONDITION OF THE SYSTEM. IF IT IS BAD DTGLS SETS NEWPT * FALSE AND RETURNS. ESTPT THEN USES DTGMOR TO GET ONE * MORE POINT AND CALLS DTGLS AGAIN, LEAVING NEWPT * FALSE. DTGLS THEN ACCUMULATES THE ADDITIONAL POINT * AND AGAIN TESTS THE CONDITION. EVENTUALLY EITHER * ENOUGH POINTS ARE ADDED TO OBTAIN ADEQUATE CONDITION * OR ELSE DTGLS APPLIES AN ARBITRARY STABILIZATION AND * SO IN EITHER CASE DTGLS FINALLY RETURNS A SOLUTION. * * ------------------------------------------------------------------ *--D replaces "?": ?TGPD, ?TGGET, ?TGMOR, ?TGLS, ?TGCM1 *--& ?TGMOR, ?TGCM1, ?TGLS, ?ROT, ?ROTG * ------------------------------------------------------------------ */ /* Note -- if nfmax is changed, mw changes and C code must change * everywhere mw is used in a dimension declarator. */ /*++ Code for .C. is ACTIVE */ /*++ Code for ~.C. is INACTIVE * double precision dsfcn, DZ(2,NP), W(mw,6) *++ END */ /* ------------------------------------------------------------------ */ #define DSFCN(x1,y1,x2,y2) ((double)(powi(((x1) + (x2))*.5 - dtgcm1.x0,2) + \ powi(((y1) + (y2))*.5 - dtgcm1.y0,2))) /* ------------------------------------------------------------------ */ dtgcm1.mfmin = 6; limit = min( MFMAX, np - 1 ); lolim = min( dtgcm1.mfmin, np - 1 ); dtgcm1.newpt = TRUE; /* SET IWORK() = 0 */ for (i = 1; i <= np; i++) { Iwork[i] = 0; } /* MAIN LOOP THRU TRIANGLES * */ it = 1; goto L_40; L_30: it += 1; L_40: if (nt < it) return; dtgget( it, trsave, triang ); /* LOOP THRU VERTICES OF A TRIANGLE */ ivert = 1; goto L_60; L_50: ivert += 1; L_60: if (ivert > 3) goto L_30; jp = Trsave[ivert + 3]; if (Iwork[jp] != 0) goto L_50; Iwork[jp] = 1; /*1001 FORMAT(39H0ESTPD.. ESTIMATING PARTIALS AT POINT ,I5/1X) * THE FORTRAN END FOLLOWS.. * WRITE(*,1001) JP * PARTIALS FOR POINT JP */ dtgcm1.x0 = X[jp]; dtgcm1.y0 = Y[jp]; dtgcm1.z0 = Z[jp]; /* SET POINTERS FOR AVAILABLE SPACE. */ for (i = 2; i <= MFMAX; i++) { dtgcm1.ke[0][i - 2] = i; } dtgcm1.avail = 1; dtgcm1.ke[0][MFMAX - 1] = 0; /* BUILD FITTING EQUATIONS. */ goto L_120; L_90: if (dtgcm1.mfit < lolim) { dtgmor( x, y, z, np, triang, w ); } /* SOLVE FITTING EQUATIONS. * */ dtgls( w, MW, dtgcm1.mfit, &dtgcm1.newpt, limit, &dz[jp - 1][0], &dz[jp - 1][1] ); /* GO FOR MORE POINTS IF NEEDED TO IMPROVE CONDITION. */ L_100: if (!dtgcm1.newpt) { dtgmor( x, y, z, np, triang, w ); dtgls( w, MW, dtgcm1.mfit, &dtgcm1.newpt, limit, &dz[jp - 1][0], &dz[jp - 1][1] ); goto L_100; } goto L_50; /* PROCESS RING OF NEIGHBORS */ L_120: for (i = 1; i <= 6; i++) { Tri[i] = Trsave[i]; } /* POP AVAIL TO L */ l = dtgcm1.avail; dtgcm1.avail = dtgcm1.ke[0][dtgcm1.avail - 1]; /* BUILD FIRST EDGE IN RING OF NEIGHBORS. */ dtgcm1.firstl = l; dtgcm1.lastl = l; dtgcm1.ke[0][l - 1] = 0; dtgcm1.ke[1][l - 1] = 0; dtgcm1.ke[2][l - 1] = it; iv1 = Add1[ivert]; dtgcm1.ke[3][l - 1] = Tri[iv1]; j1 = Tri[iv1 + 3]; iv2 = Sub1[ivert]; j2 = Tri[iv2 + 3]; Dsq[l] = DSFCN( X[j1], Y[j1], X[j2], Y[j2] ); /* BUILD EQUAS FOR FIRST TWO NEIGHBORING PTS. */ dtgcm1.mfit = 1; kp = j1; /* BUILD ROW *1000 FORMAT(25H BUILD ROW.. CENTER PT =,I5,12H, EQUA NO =,I3, * *13H NEARBY PT =,I5) * WRITE(*,1000) JP,MFIT,KP */ Jused[dtgcm1.mfit] = kp; w[3][dtgcm1.mfit - 1] = X[kp] - dtgcm1.x0; w[4][dtgcm1.mfit - 1] = Y[kp] - dtgcm1.y0; w[5][dtgcm1.mfit - 1] = Z[kp] - dtgcm1.z0; w[0][dtgcm1.mfit - 1] = SQ(w[3][dtgcm1.mfit - 1]); w[1][dtgcm1.mfit - 1] = w[3][dtgcm1.mfit - 1]*w[4][dtgcm1.mfit - 1]; w[2][dtgcm1.mfit - 1] = SQ(w[4][dtgcm1.mfit - 1]); dtgcm1.mfit = 2; kp = j2; /* BUILD ROW *1000 FORMAT(25H BUILD ROW.. CENTER PT =,I5,12H, EQUA NO =,I3, * *13H NEARBY PT =,I5) * WRITE(*,1000) JP,MFIT,KP */ Jused[dtgcm1.mfit] = kp; w[3][dtgcm1.mfit - 1] = X[kp] - dtgcm1.x0; w[4][dtgcm1.mfit - 1] = Y[kp] - dtgcm1.y0; w[5][dtgcm1.mfit - 1] = Z[kp] - dtgcm1.z0; w[0][dtgcm1.mfit - 1] = SQ(w[3][dtgcm1.mfit - 1]); w[1][dtgcm1.mfit - 1] = w[3][dtgcm1.mfit - 1]*w[4][dtgcm1.mfit - 1]; w[2][dtgcm1.mfit - 1] = SQ(w[4][dtgcm1.mfit - 1]); goto L_180; /* MOVE LEFT AND THEN RIGHT AROUND POINT JP * TO BUILD EQUAS FOR IMMEDIATELY * NEIGHBORING POINTS. * */ L_160: if (!hitbdy) goto L_90; for (i = 1; i <= 6; i++) { Tri[i] = Trsave[i]; } goto L_360; /* MOVE LEFT * LEFT MEANS COUNTERCLOCKWISE. */ L_180: jj2 = j2; /* * * **************************************** */ ipcent = ivert; hitbdy = FALSE; if (dtgcm1.mfit >= MFMAX) goto L_160; i = Sub1[ipcent]; tnow = Tri[i]; if (tnow != 0) goto L_200; /* HAVE HIT BOUNDARY */ hitbdy = TRUE; goto L_160; L_200: dtgget( tnow, tri, triang ); for (ipcent = 1; ipcent <= 3; ipcent++) { if (Tri[ipcent + 3] == jp) goto L_220; } /* IDENTIFY NEXT POINT */ L_220: i = Sub1[ipcent]; kp = Tri[i + 3]; /* RECORD INFORMATION ABOUT THE NEW EDGE IN * THE LIST STRUCTURE. * * POP AVAIL TO L */ l = dtgcm1.avail; dtgcm1.avail = dtgcm1.ke[0][dtgcm1.avail - 1]; dtgcm1.ke[1][l - 1] = dtgcm1.lastl; dtgcm1.ke[2][l - 1] = tnow; i = Add1[ipcent]; dtgcm1.ke[3][l - 1] = Tri[i]; dtgcm1.ke[0][dtgcm1.lastl - 1] = l; dtgcm1.lastl = l; jj1 = jj2; jj2 = kp; Dsq[l] = DSFCN( X[jj1], Y[jj1], X[jj2], Y[jj2] ); if (kp == j1) goto L_240; /* KP IS A NEW NEIGHBORING POINT. */ dtgcm1.mfit += 1; /* BUILD ROW *1000 FORMAT(25H BUILD ROW.. CENTER PT =,I5,12H, EQUA NO =,I3, * *13H NEARBY PT =,I5) * WRITE(*,1000) JP,MFIT,KP */ Jused[dtgcm1.mfit] = kp; w[3][dtgcm1.mfit - 1] = X[kp] - dtgcm1.x0; w[4][dtgcm1.mfit - 1] = Y[kp] - dtgcm1.y0; w[5][dtgcm1.mfit - 1] = Z[kp] - dtgcm1.z0; w[0][dtgcm1.mfit - 1] = SQ(w[3][dtgcm1.mfit - 1]); w[1][dtgcm1.mfit - 1] = w[3][dtgcm1.mfit - 1]*w[4][dtgcm1.mfit - 1]; w[2][dtgcm1.mfit - 1] = SQ(w[4][dtgcm1.mfit - 1]); dtgcm1.ke[0][l - 1] = 0; goto L_160; /* KP IS NOT A NEW POINT. WE HAVE CYCLED ALL THE WAY * AROUND TO THE INITIAL POINT J1. THIS MEANS THE * ENTIRE RING OF NEIGHBORS HAS BEEN FOUND. * */ L_240: dtgcm1.ke[0][l - 1] = dtgcm1.firstl; dtgcm1.ke[1][dtgcm1.firstl - 1] = dtgcm1.lastl; /* ABOVE LOOP IS EXITED FOR ONE OF THREE REASONS.. * (1) MFIT .EQ. MFMAX, * (2) HITBDY = .TRUE., * OR (3) COMPLETE RING OF NEIGHBORS HAS BEEN BUILT * */ goto L_160; /* ************************** * * * PROCEDURE ( MOVE RIGHT ) * * * ************************* * RIGHT MEANS CLOCKWISE. */ L_360: jj1 = j1; ipcent = ivert; L_370: if (dtgcm1.mfit >= MFMAX) goto L_90; tnow = Tri[ipcent]; /* THE FOLLOWING TEST CAUSES AN EXIT FROM THE LOOP IF THE * BOUNDARY ON THE RIGHT HAS BEEN ENCOUNTERED. * */ if (tnow == 0) goto L_90; dtgget( tnow, tri, triang ); for (ipcent = 1; ipcent <= 3; ipcent++) { if (Tri[ipcent + 3] == jp) goto L_390; } /* IDENTIFY NEXT POINT */ L_390: ip1 = Add1[ipcent]; kp = Tri[ip1 + 3]; /* RECORD INFO ABOUT THE NEW EDGE IN THE LIST STRUCTURE * * POP AVAIL TO L */ l = dtgcm1.avail; dtgcm1.avail = dtgcm1.ke[0][dtgcm1.avail - 1]; dtgcm1.ke[0][l - 1] = dtgcm1.firstl; dtgcm1.ke[1][l - 1] = 0; dtgcm1.ke[2][l - 1] = tnow; dtgcm1.ke[3][l - 1] = Tri[ip1]; dtgcm1.ke[1][dtgcm1.firstl - 1] = l; dtgcm1.firstl = l; jj2 = jj1; jj1 = kp; Dsq[l] = DSFCN( X[jj1], Y[jj1], X[jj2], Y[jj2] ); /* KP IS A NEW NEIGHBORING POINT. */ dtgcm1.mfit += 1; /* BUILD ROW *1000 FORMAT(25H BUILD ROW.. CENTER PT =,I5,12H, EQUA NO =,I3, * *13H NEARBY PT =,I5) * WRITE(*,1000) JP,MFIT,KP */ Jused[dtgcm1.mfit] = kp; w[3][dtgcm1.mfit - 1] = X[kp] - dtgcm1.x0; w[4][dtgcm1.mfit - 1] = Y[kp] - dtgcm1.y0; w[5][dtgcm1.mfit - 1] = Z[kp] - dtgcm1.z0; w[0][dtgcm1.mfit - 1] = SQ(w[3][dtgcm1.mfit - 1]); w[1][dtgcm1.mfit - 1] = w[3][dtgcm1.mfit - 1]*w[4][dtgcm1.mfit - 1]; w[2][dtgcm1.mfit - 1] = SQ(w[4][dtgcm1.mfit - 1]); goto L_370; /* ABOVE LOOP IS EXITED WHEN * (1) MFIT .EQ. MFMAX * OR (2) HIT BOUNDARY ON RIGHT * */ #undef DSFCN } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dtgmor( double x[], double y[], double z[], long np, long triang[], double w[][21]) { LOGICAL32 avok, test, tstlm1, tstlm2, tstlp1, tstlp2, used; long int i, im1, ip1, k, kp, l, limit2, lm1, lm2, lmin_, lnew, lp1, lp2, lp3, tnow, tri[6], v1, v2, _i, _r; static long int add1[3], sub1[3]; double dmin, x1, x2, y1, y2; static int _aini = 1; /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Add1 = &add1[0] - 1; double *const Dsq = &dtgcm1.dsq[0] - 1; long *const Jused = &dtgcm1.jused[0] - 1; long *const Sub1 = &sub1[0] - 1; long *const Tri = &tri[0] - 1; double *const X = &x[0] - 1; double *const Y = &y[0] - 1; double *const Z = &z[0] - 1; /* end of OFFSET VECTORS */ if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */ Add1[1] = 2; Add1[2] = 3; Add1[3] = 1; Sub1[1] = 3; Sub1[2] = 1; Sub1[3] = 2; _aini = 0; } /*>> 1995-09-26 CLL Editing for inclusion into MATH77. * * C.L.LAWSON, JPL, 1977 JAN 4 * CHANGED NOV 30, 1977 TT * THIS SUBR IS CALLED BY DTGPD AS FOLLOWS: * (1) IF NEWPT = .TRUE., * THE SUBR EXTENDS THE RING OF NEARBY POINTS AROUND * POINT (X0,Y0) UP TO A TOTAL OF MIN(MFMIN,NP-1) NEARBY POINTS * TO BE USED IN COMPUTING A LOCAL FIT TO THE SURFACE NEAR POINT * (X0,Y0). * (2) IF NEWPT = .FALSE., * THE SUBR GETS ONLY ONE MORE NEARBY POINT AND * RETURNS. * ------------------------------------------------------------------ */ /*++ Code for .C. is ACTIVE */ /*++ Code for ~.C. is INACTIVE * integer mw * parameter(mw = mfmax+5) * double precision dmin, dsfcn, W(mw,6) *++ END */ /* NAMELIST/DEBUG/AVAIL,MFMIN,MFIT,FIRSTL,LASTL, JUSED, * *X0,Y0,Z0,DSQ * ------------------------------------------------------------------ */ #define DSFCN(x1,y1,x2,y2) ((double)(powi(((x1) + (x2))*.5 - dtgcm1.x0,2) + \ powi(((y1) + (y2))*.5 - dtgcm1.y0,2))) /* ------------------------------------------------------------------ * */ if (dtgcm1.newpt) { limit2 = min( dtgcm1.mfmin, np - 1 ); } else { limit2 = dtgcm1.mfit + 1; /* WHEN DTGMOR IS CALLED WITH * NEWPT = .FALSE. MFIT MUST BE LESS THAN MFMIN * AND NP-1 SO THAT LIMIT2 CAN SAFELY BE SET * TO MFIT+1 AS ABOVE. */ } avok = TRUE; goto L_200; /* ******************************************* * * * PROCEDURE ( DELETE EDGES L,LP1, AND LP2 ) * * * ****************************************** * */ L_10: lm1 = dtgcm1.ke[1][l - 1]; lp3 = dtgcm1.ke[0][lp2 - 1]; if (lm1 != 0) dtgcm1.ke[0][lm1 - 1] = lp3; if (lp3 != 0) dtgcm1.ke[1][lp3 - 1] = lm1; if (l == dtgcm1.firstl) { dtgcm1.firstl = lp3; } else if ((lp1 == dtgcm1.firstl) || (lp2 == dtgcm1.firstl)) { dtgcm1.firstl = lp3; dtgcm1.lastl = lm1; } else if (lp2 == dtgcm1.lastl) { dtgcm1.lastl = lm1; } /* TRANSFER L, LP1, AND LP2 TO AVAILABLE SPACE LIST * */ dtgcm1.ke[0][lp2 - 1] = dtgcm1.avail; dtgcm1.avail = l; printf(" MORPTS.. DELETING THREE EDGES. TNOW=%5ld\n", tnow); goto L_200; /* * * DELETE EDGE L. ADD TWO NEW EDGES. /KP,AVOK/ */ L_50: ; /* POP AVAIL TO LNEW. SET AVOK * * * ***************************************** */ if (dtgcm1.avail == 0) { avok = FALSE; /* ERROR CONDITION * */ printf("\n SUBR DTGMOR.. WARNING.. DIMENSION OF KE(,) NOT LARGE ENOUGH TO FIND MFMIN NEARBY POINTS FOR ESTIMATING PARTIAL DERIVS.\n PROGRAM WILL CONTINUE USING FEWER NEARBY POINTS.\n"); } else { lnew = dtgcm1.avail; dtgcm1.avail = dtgcm1.ke[0][dtgcm1.avail - 1]; } if (avok) { dtgcm1.ke[0][l - 1] = lnew; dtgcm1.ke[0][lnew - 1] = lp1; dtgcm1.ke[1][lnew - 1] = l; dtgcm1.ke[2][lnew - 1] = tnow; if (lp1 != 0) dtgcm1.ke[1][lp1 - 1] = lnew; if (l == dtgcm1.lastl) dtgcm1.lastl = lnew; dtgget( tnow, tri, triang ); for (i = 1; i <= 3; i++) { if (Tri[i] == dtgcm1.ke[2][l - 1]) goto L_90; } L_90: ip1 = Add1[i]; dtgcm1.ke[3][l - 1] = Tri[ip1]; dtgcm1.ke[2][l - 1] = tnow; im1 = Sub1[i]; dtgcm1.ke[3][lnew - 1] = Tri[im1]; kp = Tri[im1 + 3]; v1 = Tri[ip1 + 3]; Dsq[l] = DSFCN( X[v1], Y[v1], X[kp], Y[kp] ); v1 = Tri[i + 3]; Dsq[lnew] = DSFCN( X[v1], Y[v1], X[kp], Y[kp] ); } /* PROCESS NEW VERTEX * * HERE THE TRIANGLE TNOW HAS ONLY EDGE L IN COMMON WITH THE * CURRENT RING AROUND POINT JP. THE VERTEX OPPOSITE SIDE L * IN TRIANGLE TNOW IS KP. * TEST TO SEE IF POINT KP HAS ALREADY BEEN USED TO * BUILD AN EQUATION. * */ used = FALSE; for (k = 1; k <= dtgcm1.mfit; k++) { if (Jused[k] == kp) { used = TRUE; goto L_170; } } L_170: if (used) goto L_200; dtgcm1.mfit += 1; /* BUILD ROW *1000 FORMAT(' BUILD ROW IN DTGMOR.. ',5X,'' EQUA NO =',I3, * * ' NEARBY PT =',I5) * WRITE(*,1000) MFIT,KP */ Jused[dtgcm1.mfit] = kp; w[3][dtgcm1.mfit - 1] = X[kp] - dtgcm1.x0; w[4][dtgcm1.mfit - 1] = Y[kp] - dtgcm1.y0; w[5][dtgcm1.mfit - 1] = Z[kp] - dtgcm1.z0; w[0][dtgcm1.mfit - 1] = SQ(w[3][dtgcm1.mfit - 1]); w[1][dtgcm1.mfit - 1] = w[3][dtgcm1.mfit - 1]*w[4][dtgcm1.mfit - 1]; w[2][dtgcm1.mfit - 1] = SQ(w[4][dtgcm1.mfit - 1]); L_200: if ((dtgcm1.mfit >= limit2) || (!avok)) return; /* WRITE(*,DEBUG) *1004 FORMAT(/' KE(,)='/(7X,5I5)) * WRITE(*,1004) (I,(KE(I,J), J=1,4), I=1,16) * * FIND EDGE WHOSE MIDPOINT IS CLOSEST TO POINT JP. * SKIP BOUNDARY EDGES. */ l = dtgcm1.firstl; lmin_ = 0; L_220: if (dtgcm1.ke[3][l - 1] != 0) { if (lmin_ == 0) { test = TRUE; } else { test = Dsq[l] < dmin; } if (test) { dmin = Dsq[l]; lmin_ = l; } } if (l == dtgcm1.lastl) goto L_230; l = dtgcm1.ke[0][l - 1]; goto L_220; L_230: l = lmin_; /* EDGE L IS THE NEAREST TO POINT JP. * INVESTIGATE TRIANGLE ON OPPOSIDE SIDE OF EDGE L. * */ tnow = dtgcm1.ke[3][l - 1]; lp1 = dtgcm1.ke[0][l - 1]; tstlp1 = lp1 != 0; if (tstlp1) tstlp1 = dtgcm1.ke[3][lp1 - 1] == tnow; lm1 = dtgcm1.ke[1][l - 1]; tstlm1 = lm1 != 0; if (tstlm1) tstlm1 = dtgcm1.ke[3][lm1 - 1] == tnow; if (!tstlp1) goto L_270; lp2 = dtgcm1.ke[0][lp1 - 1]; tstlp2 = lp2 != 0; if (tstlp2) tstlp2 = dtgcm1.ke[3][lp2 - 1] == tnow; if (!tstlp2) goto L_250; goto L_10; L_250: if (!tstlm1) goto L_300; lp2 = lp1; lp1 = l; l = lm1; goto L_10; L_270: if (!tstlm1) goto L_50; lm2 = dtgcm1.ke[1][lm1 - 1]; tstlm2 = lm2 != 0; if (tstlm2) tstlm2 = dtgcm1.ke[3][lm2 - 1] == tnow; if (!tstlm2) goto L_280; lp2 = l; lp1 = lm1; l = lm2; goto L_10; L_280: lp1 = l; l = lm1; /* END OF TOP LEVEL CODE. PROCEDURES FOLLOW.. * * **************************************************** * * * PROCEDURE ( DELETE EDGES L AND LP1. ADD NEW EDGE) * * * *************************************************** * */ L_300: lp2 = dtgcm1.ke[0][lp1 - 1]; dtgcm1.ke[0][l - 1] = lp2; if (lp2 != 0) dtgcm1.ke[1][lp2 - 1] = l; if (lp1 == dtgcm1.firstl) { dtgcm1.firstl = lp2; } else if (lp1 == dtgcm1.lastl) { dtgcm1.lastl = l; } dtgcm1.ke[0][lp1 - 1] = dtgcm1.avail; dtgcm1.avail = lp1; dtgget( tnow, tri, triang ); for (i = 1; i <= 3; i++) { if (Tri[i] == dtgcm1.ke[2][l - 1]) goto L_320; } L_320: dtgcm1.ke[2][l - 1] = tnow; ip1 = Add1[i]; dtgcm1.ke[3][l - 1] = Tri[ip1]; if (dtgcm1.ke[3][l - 1] != 0) { im1 = Sub1[i]; v1 = Tri[ip1 + 3]; v2 = Tri[im1 + 3]; Dsq[l] = DSFCN( X[v1], Y[v1], X[v2], Y[v2] ); } /*1003 FORMAT(60X,'MORPTS.. DELETING TWO EDGES, ADDING ONE.', * * ' TNOW=', I5) * WRITE(*,1003) TNOW */ goto L_200; #undef DSFCN } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dtgls( double w[][21], long mw, long mfit, LOGICAL32 *newpt, long limit, double *dz1, double *dz2) { long int i, irow, j, jj, kp, len; double avesq, c1, c2, cosine, dmin, sine, stab, sum; static double ave, wt, wt1, wt2; /*>> 1995-09-26 CLL Editing for inclusion into MATH77. * * WEIGHT, SCALE, STABILIZE, AND SOLVE A LEAST SQRS PROBLEM. * * C.L.LAWSON, JPL, 1976 FEB 27 * CHANGED NOV 28, 1977 TT * EDITED CODE. NO CHANGE IN ALGORITHM. CLL. 1979 MAR 5. * * ON INPUT, WITH NEWPT TRUE, COLUMNS OF W(,) CONTAIN TRANSLATED * VALUES OF X**2, X*Y, Y**2, X, Y, Z. * * THE NUMBER, STAB, IS USED AS A MARQUARDT STABILIZATION * VALUE TO DAMP DOWN THE VALUES OF THE THREE SECOND PARTIAL * DERIVATIVES WHEN THE AVAILABLE DATA DOES NOT PRODUCE A * REASONABLY WELL-CONDITIONED LEAST SQUARES PROBLEM. * A LARGER VALUE OF STAB DAMPS THE SECOND PARTIALS MORE. * * THE OUTPUT IS DZ1 AND DZ2. * ------------------------------------------------------------------ */ /*++ Code for .C. is ACTIVE */ /*++ Code for ~.C. is INACTIVE * double precision SINE, STAB, SUM, W(MW,6) *++ END */ /* ------------------------------------------------------------------ */ stab = 1.e0; if (!*newpt) goto L_200; sum = 0.0e0; kp = mfit; for (i = 1; i <= kp; i++) { sum += w[0][i - 1] + w[2][i - 1]; } /* KP will be converted to floating point. */ avesq = sum/kp; ave = sqrt( avesq ); /* WE HAVE EXPERIMENTED WITH A VARIETY OF * WEIGHTING FORMULAS AND FOUND NONE SIG- * NIFICANTLY BETTER THAN THE UNIFORM * WEIGHTING.. WT = 1. */ wt = 1.0e0; wt1 = wt/ave; wt2 = wt/avesq; irow = 1; goto L_130; /* SCALE ROW IROW */ L_120: w[0][irow - 1] *= wt2; w[1][irow - 1] *= wt2; w[2][irow - 1] *= wt2; w[3][irow - 1] *= wt1; w[4][irow - 1] *= wt1; w[5][irow - 1] *= wt; irow += 1; L_130: if (kp >= irow) goto L_120; /* GIVENS TRIANGULARIZATION */ for (j = 1; j <= 6; j++) { if (j + 1 <= kp) { len = 6 - j; for (i = j + 1; i <= kp; i++) { if (w[j - 1][i - 1] != 0.0e0) { drotg( &w[j - 1][j - 1], &w[j - 1][i - 1], &cosine, &sine ); w[j - 1][i - 1] = 0.0e0; if (len > 0) drot( len, &w[j][j - 1], mw, &w[j][i - 1], mw, cosine, sine ); } } } } goto L_230; L_200: irow = mfit; /* SCALE NEW ROW */ w[0][irow - 1] *= wt2; w[1][irow - 1] *= wt2; w[2][irow - 1] *= wt2; w[3][irow - 1] *= wt1; w[4][irow - 1] *= wt1; w[5][irow - 1] *= wt; /* THE NEW ROW WILL BE ACCUMULATED INTO THE TRIANGLE * USING GIVENS ROTATIONS. */ for (j = 1; j <= 6; j++) { len = 6 - j; if (w[j - 1][irow - 1] != 0.0e0) { drotg( &w[j - 1][j - 1], &w[j - 1][irow - 1], &cosine, &sine ); w[j - 1][irow - 1] = 0.0e0; if (len > 0) drot( len, &w[j][j - 1], mw, &w[j][irow - 1], mw, cosine, sine ); } } L_230: if (kp >= 5) { dmin = fmin( fmin( fmin( fabs( w[0][0] ), fabs( w[1][1] ) ), fmin( fabs( w[2][2] ), fabs( w[3][3] ) ) ), fabs( w[4][4] ) ); } else { dmin = 0.0e0; } /* TEST CONDITION OF SYSTEM. */ if (dmin >= 0.01e0) goto L_300; /* SYSTEM IS ILL-CONDITIONED. * */ if (mfit < limit) { /* RETURN TO ADD ONE MORE POINT TO THE FITTING PROBLEM * IN HOPES THAT THIS WILL IMPROVE THE CONDITION. * */ *newpt = FALSE; return; } /* CANNOT ADD ANY MORE POINTS. STABILIZE BY DAMPING * THE SECOND PARTIAL DERIVATIVES. * */ jj = 1; goto L_280; L_250: ; /* THE NEW ROW WILL BE ACCUMULATED INTO THE TRIANGLE * USING GIVENS ROTATIONS. */ for (j = 1; j <= 6; j++) { len = 6 - j; if (w[j - 1][irow - 1] != 0.0e0) { drotg( &w[j - 1][j - 1], &w[j - 1][irow - 1], &cosine, &sine ); w[j - 1][irow - 1] = 0.0e0; if (len > 0) drot( len, &w[j][j - 1], mw, &w[j][irow - 1], mw, cosine, sine ); } } jj += 1; L_280: if (jj > 3) goto L_300; irow = mfit + jj; for (j = 1; j <= 6; j++) { w[j - 1][irow - 1] = 0.0e0; } w[jj - 1][irow - 1] = stab; goto L_250; /* SOLVE FOR DZ1 AND DZ2 * */ L_300: c2 = w[5][4]/w[4][4]; c1 = (w[5][3] - c2*w[4][3])/w[3][3]; *dz1 = c1/ave; *dz2 = c2/ave; *newpt = TRUE; return; } /* end of function */