/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:32:03 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dsquad.h" #include /* PARAMETER translations */ #define HALF 0.50e0 #define ZERO 0.0e0 /* end of PARAMETER translations */ double /*FUNCTION*/ dsquad( long k, long n, double t[], double bcoef[], double x1, double x2) { long int i, jf, left, m, mf, mode; double a, aa, b, bb, bma, bpa, c1, dsquad_v, q, sum[5], ta, tb; static long il1 = 1; static long il2 = 1; static double gpts[9]={0.5773502691896257645091487805019574556476e0, 0.2386191860831969086305017216807119354186e0,0.6612093864662645136613995950199053470064e0, 0.9324695142031520278123015544939946091348e0,0.1488743389816312108848260011297199846176e0, 0.4333953941292471907992659431657841622001e0,0.6794095682990244062343273651148735757693e0, 0.8650633666889845107320966884234930485275e0,0.9739065285171717200779640120844520534283e0}; static double gwts[9]={1.0e0,0.4679139345726910473898703439895509948120e0, 0.3607615730481386075698335138377161116612e0,0.1713244923791703450402961421727328935273e0, 0.2955242247147528701738929946513383294210e0,0.2692667193099963550912269215694693528586e0, 0.2190863625159820439955349342281631924634e0,0.1494513491505805931457763396576973323908e0, 0.06667134430868813759356880989333179285686e0}; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Bcoef = &bcoef[0] - 1; double *const Gpts = &gpts[0] - 1; double *const Gwts = &gwts[0] - 1; double *const Sum = &sum[0] - 1; double *const T = &t[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. *>> 1996-03-30 DSQUAD Krogh Added external statement. *>> 1994-10-20 DSQUAD Krogh Changes to use M77CON *>> 1992-11-12 DSQUAD C. L. Lawson, JPL Saving IL1 & IL2. *>> 1992-10-27 C. L. Lawson, JPL *>> 1988-04-01 C. L. Lawson, JPL * * This subprogram computes the integral over [X1,X2] of a K-th order * B-spline using the B-representation (T,BCOEF,N,K). Orders * K as high as 20 are permitted. Permits X1 .le. X2 or X1 gt. X2. * The "proper interpolation interval" is from T(K) to T(N+1). * For most reliable results, X1 and X2 should lie in this * interval, however this subr will give a result even if this * is not the case by use of extrapolation. * This subr applies a 2, 6, or 10 point Gauss formula on * subintervals of [X1,X2] which are formed by the included * (distinct) knots. */ /* DESCRIPTION OF ARGUMENTS * INPUT * K - ORDER OF B-SPLINE, 2 .le. K .le. 20 * Note that the polynomial degree of the segments of * the spline are one less than the order. * N - LENGTH OF COEFFICIENT ARRAY * T() - KNOT ARRAY OF LENGTH N+K * BCOEF() - B-SPLINE COEFFICIENT ARRAY OF LENGTH N * X1,X2 - END POINTS OF QUADRATURE INTERVAL. * Permit X1 .le. X2 or X1 .gt. X2. * OUTPUT * DSQUAD - INTEGRAL OF THE B-SPLINE OVER [X1,X2]. * * ERROR CONDITIONS: None * ------------------------------------------------------------------ * Original code was called BQUAD or BSQUAD. Written by D. E. Amos, * Sandia, June, 1979. Documented in SAND79-1825. * Uses the (T, BCOEF, N, K) representation of a B-spline as * presented in A PRACTICAL GUIDE TO SPLINES by Carl De Boor, * Springer-Verlag, 1978. * Current version by C. L. Lawson, JPL, March 1988. * ------------------------------------------------------------------ *--D replaces "?": ?SQUAD, ?SFIND, ?SVAL * ------------------------------------------------------------------ * GPTS(), GWTS() Index 1 selects 2-point Gaussian formula. * Indices 2-4 select 6-point Gaussian formula. * Indices 5-9 select 10-point Gaussian formula. * ------------------------------------------------------------------ */ /* DATA GPTS / 5.77350269189626D-01, 2.38619186083197D-01, * 1 6.61209386466265D-01, 9.32469514203152D-01, 1.48874338981631D-01, * 2 4.33395394129247D-01, 6.79409568299024D-01, 8.65063366688985D-01, * 3 9.73906528517172D-01/ * * DATA GWTS / 1.00000000000000D+00, 4.67913934572691D-01, * 1 3.60761573048139D-01, 1.71324492379170D-01, 2.95524224714753D-01, * 2 2.69266719309996D-01, 2.19086362515982D-01, 1.49451349150581D-01, * 3 6.66713443086880D-02/ * * Gaussian quadrature abcissae and weights to 40 decimal digits. */ /* ------------------------------------------------------------------ * */ aa = fmin( x1, x2 ); bb = fmax( x1, x2 ); if (aa == bb) { q = ZERO; goto L_90; } /* Selection of Gaussian quadrature formula: * KORDER .le. 4 => 2 point formula * 5 .le. KORDER .le. 12 => 6 point formula * KORDER .ge. 13 => 10 point formula * */ if (k <= 4) { jf = 0; mf = 1; } else if (k <= 12) { jf = 1; mf = 3; } else { jf = 4; mf = 5; } for (i = 1; i <= mf; i++) { Sum[i] = ZERO; } dsfind( t, k, n + 1, aa, &il1, &mode ); dsfind( t, k, n + 1, bb, &il2, &mode ); for (left = il1; left <= il2; left++) { ta = T[left]; tb = T[left + 1]; if (ta != tb) { if (left == k) { a = aa; } else { a = fmax( aa, ta ); } if (left == n) { b = bb; } else { b = fmin( bb, tb ); } bma = HALF*(b - a); bpa = HALF*(b + a); for (m = 1; m <= mf; m++) { c1 = bma*Gpts[jf + m]; Sum[m] += bma*(dsval( k, n, t, bcoef, bpa + c1, 0 ) + dsval( k, n, t, bcoef, bpa - c1, 0 )); } } } q = ZERO; for (m = 1; m <= mf; m++) { q += Gwts[jf + m]*Sum[m]; } if (x1 > x2) q = -q; L_90: ; dsquad_v = q; return( dsquad_v ); } /* end of function */