/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:57 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "sc2fit.h" #include /* PARAMETER translations */ #define NBAND 4 #define NBAND1 (NBAND + 1) #define ONE 1.0e0 #define ZERO 0.0e0 /* end of PARAMETER translations */ void /*FUNCTION*/ sc2fit( float xi[], float yi[], float sdi[], long nxy, float b[], long nb, float *w, long nw, float yknot[], float ypknot[], float *sigfac, long *ierr1) { #define W(I_,J_) (*(w+(I_)*(nw)+(J_))) LOGICAL32 newseg, usewt1; long int i, ierr2, ierr3, irnow, iseg, j, jpoint, jtprev, k, ksize, n1, nparam; float dof, p[4], rnorm, sdijp, wt, wt1; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const B = &b[0] - 1; float *const P = &p[0] - 1; float *const Sdi = &sdi[0] - 1; float *const Xi = &xi[0] - 1; float *const Yi = &yi[0] - 1; float *const Yknot = &yknot[0] - 1; float *const Ypknot = &ypknot[0] - 1; /* end of OFFSET VECTORS */ /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 2000-12-01 SC2FIT Krogh Dim. SDI(*) instead NXY. *>> 1995-11-21 SC2FIT Krogh Converted from SFTRAN to Fortran 77. *>> 1994-10-19 SC2FIT Krogh Changes to use M77CON *>> 1994-01-31 SC2FIT CLL Added test for SDI(i) .le. 0 when SDI(1) > 0. *>> 1990-01-23 CLL Deleted ref to unused variable NX in call to IERM1 *>> 1989-10-20 CLL *>> 1987-10-22 SC2FIT Lawson Initial code. * Least squares fit to discrete data by a C-2 cubic spline. * ------------------------------------------------------------------ * Algorithm and program designed by C.L.Lawson and R.J.Hanson. * The general approach but not the complete code is given in * 'SOLVING LEAST SQUARES PROBLEMS', by Lawson and Hanson, * publ by Prentice-Hall, 1974. * Programming and later changes and corrections by Lawson,Hanson, * T.Lang, and D.Campbell, Sept 1968, Nov 1969, and Aug 1970. * Modified 1968 Sept 17 to provide C-2 continuity. * 1974 5/21, C.L.Lawson, Fixed bug that caused ISEG to get too big. * Also changed to exit immidiately if B() array is not * strictly increasing. * 1984 July 10. Modified for improved portability and to conform * to Fortran 77. C. L. Lawson, JPL. * Added calls to the error message subrs. * 7/23/87 CLL. Added the IERR1 argument. * 1989-10-20 CLL Changed code so there are no RETURN statements in * Sftran procedures. Previous code had such a RETURN that led to * a warning diagnostic from a Cray compiler due to unreachable * CONTINUE statements. Also introduced "c--" lines for CHGTYP. * ------------------------------------------------------------------ * SUBROUTINE ARGUMENTS * * (XI(i),i=1,NXY) [in] Abcissas of data to be fitted. Require * this data be ordered so X(i) .le. X(i+1). * * (YI(i),i=1,NXY) [in] Ordinates of data to be fitted. * * (SDI(i),i=1,NXY) [in] User may use this array to assign an * a priori standard deviation of error to each * YI(i) value. The weighted fitting algorithm will take * account of these. Optionally the user may set SDI(1) to * a negative value. Then this subr will use ABS(SDI(1)) as * the standard deviation for each YI(i) value. In this case * the SDI() array can be dimensioned SDI(1). * If SDI(1) = 0., the subr issues an error message and returns. * * NXY [in] No. of data pairs, (XI(i), YI(i)), and no. of elts * in SDI() if SDI(1) is positive. Require NXY .ge. 4. * * (B(j),j=1,NB) [in] Breakpoints for the spline function, * including endpoints. These breakpoints must be * strictly increasing: B(j) .lt. B(j+1). * It is required that all abcissas, XI(i), lie in the * closed interval, [B(1), B(NB)]. * * NB [in] No. of breakpoints, including endpoints. * The no. of parameters in the least squares problem * will be NB + 2. * To have a nonsingular problem one must have * NXY .ge. NB + 2, and the distribution of the breakpoints * must not be too skewed relative to the data abcissas. * If singularity is detected, an error message will be * issued by the subr that solves the band matrix. * * W() [scratch] Work space dimensioned W(NW,5). * * NW [in] First dimension of W(). Must satisfy NW .ge. NB + 4 * Let KMAX denote the max no. of data abcissas, XI(i), * in any one breakpoint interval, i.e. between B(j) and * B(j+1) for some j. The subr will be more efficient * if NW is at least NB + 3 + KMAX. * * (YKNOT(k) and YPKNOT(k),k=1,NB) [out] The subr will return * values defining the fitted C2 spline curve in these arrays. * These values and first derivatives of the fitted curve at the * knot abcissae. YKNOT(j) = f(B(j)) and * YPKNOT(j) = fprime(B(j)) for j = 1,...,NB. * The user can then evaluate the fitted curve at any point by * Hermite interpolation. See subrs DHINT or SHINT. * * SIGFAC [out] The subr sets SIGFAC to RNORM / sqrt(DOF) where * RNORM = sqrt( sum over i of [( (yfit(i) - YI(i))/SDI(i))**2]) * and DOF = max(1, NXY - (NB+2)) * * IERR1 [out] Error status indicator. Note that IERR2 comes from * SBACC and IERR3 comes from SBSOL. * * = 0 means no errors detected. * = 100 means NB .lt. 2 .or. NXY .lt. NB+2 * = 200 means B(I) .ge. B(I+1) * = 300 means NW .lt. NB+4 * = 400 means XI(I-1) .gt. XI(I) * = 500 means B(1) .gt. XI(1) .or. B(NB) .lt. XI(NXY) * = 600 means Need larger dimension NW. * = 700 + IERR2 means IERR2 .ne. 0 * = 800 + IERR2 means IERR2 .ne. 0 * = 900 + IERR2 means IERR2 .ne. 0 * = 1000 + IERR3 means IERR3 .ne. 0 due to singularity * detected in SBSOL. */ /* = 1100 means SDI(1) = zero. * = 1200 means SDI(1) > zero and SDI(i) .le. zero for some i. * ------------------------------------------------------------------ * Important internal variables. * * ISEG Index of current spline segment, starting with 1 for * the first segment. * Also tells the band matrix subroutine the column index * of the least squares matrix with which the first col * of the new block of data in G() is to be associated. * KSIZE Size of current block. * JPOINT Current data pointer. * ------------------------------------------------------------------ *--S replaces "?": ?C2FIT, ?BACC, ?C2BAS, ?ERM1, ?ERV1, ?BSOL, ?TRC2C * Both versions use ERMSG, IERM1, IERV1 * Lower level subrs needed: (D/S)HTCC, (D/S)NRM2, ERFIN * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ *ierr1 = 0; /* EXIT IMMEDIATELY IF NB .lt. 2 OR NXY .LT NB+2 OR IF THE * BREAKPOINTS ARE NOT STRICTLY INCREASING. * */ nparam = nb + 2; if (nb < 2 || nxy < nparam) { *ierr1 = 100; ierm1( "SC2FIT", *ierr1, 0, "Require NB .ge. 2 and NXY .ge. NB+2" , "NB", nb, ',' ); ierv1( "NXY", nxy, '.' ); goto L_300; } n1 = nb - 1; for (i = 1; i <= n1; i++) { if (B[i] >= B[i + 1]) { *ierr1 = 200; ierm1( "SC2FIT", *ierr1, 0, "Require knots, B(I), to be strictly increasing." , "I", i, ',' ); serv1( "B(I)", B[i], ',' ); serv1( "B(I+1)", B[i + 1], '.' ); goto L_300; } } /* Require NW .ge. NB+4 * */ if (nw < nb + 4) { *ierr1 = 300; ierm1( "SC2FIT", *ierr1, 0, "Require NW .ge. NB+4", "NW", nw, ',' ); ierv1( "NB", nb, '.' ); goto L_300; } /* ------------------------------------------------------------------ * TEST SDI(1) */ if (Sdi[1] < ZERO) { wt1 = -ONE/Sdi[1]; usewt1 = TRUE; } else if (Sdi[1] > ZERO) { usewt1 = FALSE; } else { *ierr1 = 1100; ermsg( "SC2FIT", *ierr1, 0, "Require SD(1) .ne. Zero", '.' ); return; } /* Test ordering of XI() array. * */ for (i = 2; i <= nxy; i++) { if (Xi[i - 1] > Xi[i]) { *ierr1 = 400; ierm1( "SC2FIT", *ierr1, 0, "Require abcissas, X(I), to be nondecreasing." , "I", i, ',' ); serv1( "X(I-1)", Xi[i - 1], ',' ); serv1( "X(I)", Xi[i], '.' ); goto L_300; } } /* TEST THE FIRST AND LAST BREAKPOINT * FOR BRACKETING THE DATA ABCISSAS. * */ if (B[1] > Xi[1] || B[nb] < Xi[nxy]) { *ierr1 = 500; serm1( "SC2FIT", *ierr1, 0, "Require B(1) .LE. X(1) and B(NB) .ge. XI(NXY)" , "B(1)", B[1], ',' ); serv1( "X(1)", Xi[1], ',' ); serv1( "B(NB)", B[nb], ',' ); serv1( "X(NXY)", Xi[nxy], '.' ); goto L_300; } /* BEGIN LOOP TO FORM EQUATIONS FOR C2 LEAST SQUARES FIT. * */ irnow = 1; k = 1; ksize = 0; newseg = TRUE; iseg = 1; for (jpoint = 1; jpoint <= nxy; jpoint++) { if (k > nw) { sbacc( w, nw, NBAND, &irnow, ksize, iseg, &jtprev, &ierr2 ); if (ierr2 != 0) { *ierr1 = 700 + ierr2; goto L_200; } if (irnow > nw) { *ierr1 = 600; ierm1( "SC2FIT", *ierr1, 0, "Need larger dimension NW." , "NW", nw, '.' ); goto L_300; } k = irnow; ksize = 0; } /* DO WHILE( XI(JPOINT) .gt. B(ISEG+1) ) */ L_80: if (Xi[jpoint] > B[iseg + 1]) { sbacc( w, nw, NBAND, &irnow, ksize, iseg, &jtprev, &ierr2 ); if (ierr2 != 0) { *ierr1 = 800 + ierr2; goto L_200; } ksize = 0; k = irnow; iseg += 1; newseg = TRUE; goto L_80; } /* END WHILE * Build one equation */ sc2bas( Xi[jpoint], 0, iseg, &newseg, b, nb, p ); if (usewt1) { wt = wt1; } else { sdijp = Sdi[jpoint]; if (sdijp > ZERO) { wt = ONE/sdijp; } else { *ierr1 = 1200; ermsg( "SC2FIT", *ierr1, 0, "With SD(1) > 0 require all SD(I) > 0." , ',' ); serv1( "SD(1)", Sdi[1], '.' ); ierv1( "I", jpoint, ',' ); serv1( "SD(I)", sdijp, '.' ); return; } } for (j = 1; j <= 4; j++) { W(j - 1,k - 1) = P[j]*wt; } W(4,k - 1) = Yi[jpoint]*wt; /* End of build one equation */ k += 1; ksize += 1; } sbacc( w, nw, NBAND, &irnow, ksize, iseg, &jtprev, &ierr2 ); if (ierr2 != 0) { *ierr1 = 900 + ierr2; goto L_200; } /* ALL DATA POINTS HAVE BEEN PROCESSED. CALL FOR SOLUTION. * */ sbsol( 1, w, nw, NBAND, irnow, jtprev, &W(NBAND1 - 1,0), nparam, &rnorm, &ierr3 ); if (ierr3 != 0) { *ierr1 = 1000 + ierr2; ermsg( "SC2FIT", *ierr1, 0, "Singularity noted in SBSOL.", '.' ); goto L_300; } dof = max( 1, nxy - nparam ); *sigfac = rnorm/sqrtf( dof ); /* TRANSFORM PARAMETERS TO Y,YPRIME BASIS * */ strc2c( b, nb, &W(4,0), yknot, ypknot ); return; /* Error in _BACC */ L_200: ermsg( "SC2FIT", *ierr1, 0, "Error detected in subroutine SBACC" , '.' ); /* Set YKNOT() & YPKNOT() to zero */ L_300: for (i = 1; i <= nb; i++) { Yknot[i] = ZERO; Ypknot[i] = ZERO; } return; #undef W } /* end of function */