/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:49 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "swcomp.h" #include /* PARAMETER translations */ #define NMAX 10 #define NMAXP1 (NMAX + 1) /* end of PARAMETER translations */ void /*FUNCTION*/ swatan( long n5, float u5[], float y5[]) { long int i, ns; float s1[NMAXP1]; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const S1 = &s1[0] - 1; float *const U5 = &u5[0] - 1; float *const Y5 = &y5[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. *>> 2001-06-18 SWCOMP Krogh Replaced ". and." with " .and." *>> 1996-04-15 SWCOMP Krogh Changed single common block name to SWCOM. *>> 1995-11-21 SWCOMP Krogh Removed multiple entries for C conversion. *>> 1994-11-11 SWCOMP Krogh Declared all vars. *>> 1994-10-31 SWCOMP Krogh Changes to use M77CON *>> 1987-12-07 SWCOMP Lawson Initial code. *--S replaces "?": ?WCOMP,?WATAN,?WASIN,?WACOS,?WATN2,?WSUM,?WDIF, *--& ?WSQRT,?WEXP,?WSIN,?WCOS,?WTAN,?WSINH,?WCOSH,?WTANH,?WSET,?WSUM1, *--& ?WDIF1,?WPRO1,?WQUO1,?WLOG,?WPWRI,?WCHN,?WPRO,?WQUO,?PASCL, *--& ?WRCHN,?WACSI,?WCOM * * File SWCOMP [- Wengert derivative COMPutation] * The file SWCOMP contains a set of program units to perform * computation on (N+1)-tuples representing the value of a quantity * and its first N derivatives with respect to a single variable. * In the comments we use T as the name of the ultimate independent * variable. * The parameter NMAX must be set in program units * SWATAN, SWSUM, and SWPRO * to establish the largest order of derivative the package * will be able to handle. * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * C. L. Lawson, JPL, 1971. * References: * 1. Wengert, R. E., A simple automatic derivative evaluation * program, Comm. ACM, 1, Aug 1964, 463-464. * 2. C. L. Lawson, Computing Derivatives using W-arithmetic and * U-arithmetic., JPL Appl Math TM 289, Sept 1971. * Revised by CLL for Fortran 77 in Jan 1987. Deleted subr WSTART. * Put a first time flag in WPRO. Now user does not need to make * any initialization call. * 9/18/87, CLL. Added new entry names SWASIN & SWACOS. * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * NMAX Setting NMAX in the parameter statement determines * the highest order derivative this program unit will be able * to handle. The following arrays have dimensions depending * on NMAX: S1(), S2(), S3(), S4(), C() * * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * ERROR HANDLING * * All detected error conditions are essentially fatal for the * requested operation. We do not stop; rather, we issue an error * message and return. By use of the the MATH77 library subroutine, * ERMSET, the user can change the error processing action to cause a * STOP following the error message. * * Error Called * No. name Explanation * * 1 SWASIN Infinite deriv when arg = -1 or +1 * 1 SWACOS Infinite deriv when arg = -1 or +1 * 2 SWSQRT Infinite deriv when arg = 0 * 3 SWQUO1 Zero divisor * 4 SWPWRI U**M is infinite when U = 0 and M < 0 * 5 SWPRO Require dimension NMAX .ge. NDERIV * 6 SWQUO Require dimension NMAX .ge. NDERIV * 7 SWQUO Zero divisor * 8 SPASCL Require N .ge. 0 * 9 SWRCHN Require U(2) .ne. 0. * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * PROGRAM UNITS, entry NAMES, and calls * * All of these subroutine and entry names may be called directly by * a user, however, it is expected that SWCHN, SWRCHN, and SPASCL * would only be used to augment the package with new functions. * * Subroutine names: Calls to: * * subroutine SWATAN (N5,U5,Y5) SWPWRI,SWCHN * subroutine SWASIN (N16,U16,Y16) SWACSI * subroutine SWACOS (N16,U16,Y16) SWACSI * subroutine SWACSI (N16,U16,Y16,ACOSIN) SWSQRT,SWPWRI,SWPRO1,SWCHN * subroutine SWATN2 (N9,Y9,X9,A9) SWSUM,SWDIF,SWPRO,SWQUO,SWCHN * * subroutine SWSUM (NDERIV,U,V,W) None * subroutine SWDIF (N2,U2,V2,W2) None * subroutine SWSQRT (N15,U15,Y15) SWCHN * subroutine SWEXP (N4,U4,Y4) SWCHN * subroutine SWSIN (N7,U7,Y7) SWCHN * subroutine SWCOS (N7,U7,Y7) SWCHN * subroutine SWTAN (N7,U7,Y7) SWCHN * subroutine SWSINH (N7,U7,Y7) SWCHN * subroutine SWCOSH (N7,U7,Y7) SWCHN * subroutine SWTANH (N7,U7,Y7) SWCHN, SWPRO, SWQUO * subroutine SWSET (N10,UVAL,UDER,U10) None * subroutine SWSUM1 (N11,C11,U11,Y11) None * subroutine SWDIF1 (N12,C12,U12,Y12) None * subroutine SWPRO1 (N13,C13,U13,Y13) None * subroutine SWQUO1 (N14,C14,U14,Y14) SWQUO * * subroutine SWLOG (NDERIV,U,Y) SWCHN * subroutine SWPWRI (NDERIV,MPWR,U,Y) SWCHN * * subroutine SWCHN (NDERIV,U,F) SWPRO * * subroutine SWPRO (NDERIV,U,V,W) SPASCL */ /* subroutine SWQUO (NDERIV,U,V,W) SPASCL * * subroutine SPASCL (N ,C) None * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Dependencies between program units * * The progran units are SWATAN,SWSUM,SWLOG,SWCHN,SWPRO,SPASCL * SWATAN has calls into SWSUM,SWLOG,SWCHN,SWPRO * SWSUM has calls into SWLOG,SWCHN,SWPRO * SWLOG has calls into SWCHN * SWCHN has calls into SWPRO * SWPRO has calls into SPASCL * SPASCL has no calls * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * * COMPUTE.. Y5=ATAN(U5) and derivs w.r.t. T * */ Y5[1] = atanf( U5[1] ); if (n5 == 0) return; /* deriv of Y w.r.t. U IS 1./(1. + U**2) * START BY CONSTRUCTING 1.+U**2 AND ITS derivs IN S1(). * */ S1[1] = 1.0e0 + SQ(U5[1]); ns = n5 - 1; if (ns == 0) goto L_52; S1[2] = U5[1] + U5[1]; if (ns == 1) goto L_52; S1[3] = 2.0e0; if (ns == 2) goto L_52; for (i = 3; i <= ns; i++) { S1[i + 1] = 0.0e0; } /* NOW S1() = 1.+U**2 and derivs w.r.t. U . * COMPUTE 1./S1 and derivs w.r.t. U. STORE RESULT STARTING IN Y(2) * */ L_52: swpwri( ns, -1, s1, &Y5[2] ); /* Convert Y from derivs w.r.t. U TO derivs w.r.t. T */ swchn( n5, u5, y5 ); return; } /* end of function */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ void /*FUNCTION*/ swasin( long n16, float u16[], float y16[]) { /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U16 = &u16[0] - 1; float *const Y16 = &y16[0] - 1; /* end of OFFSET VECTORS */ /* COMPUTES.. Y = ASIN(U) and derivs w.r.t. T * */ Y16[1] = asinf( U16[1] ); if (n16 != 0) swacsi( n16, u16, y16, FALSE ); return; } /* end of function */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ void /*FUNCTION*/ swacos( long n16, float u16[], float y16[]) { /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U16 = &u16[0] - 1; float *const Y16 = &y16[0] - 1; /* end of OFFSET VECTORS */ /* COMPUTES.. Y = ACOS(U) and derivs w.r.t. T * */ Y16[1] = acosf( U16[1] ); if (n16 != 0) swacsi( n16, u16, y16, TRUE ); return; } /* end of function */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ void /*FUNCTION*/ swacsi( long n16, float u16[], float y16[], LOGICAL32 acosin) { long int i, ns; float s1[NMAXP1], s2[NMAXP1]; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const S1 = &s1[0] - 1; float *const S2 = &s2[0] - 1; float *const U16 = &u16[0] - 1; float *const Y16 = &y16[0] - 1; /* end of OFFSET VECTORS */ /* Common code for SWACOS and SWASIN -- ACOSIN - .true. for SWACOS */ /* Deriv of Y w.r.t. U is (+,-)1.0/sqrt(1.0 - U**2) * Start by constructing 1.0 - U**2 and its derivs in S1(). * */ S1[1] = 1.0e0 - SQ(U16[1]); if (S1[1] == 0.0e0) { /* Error condition. */ if (acosin) { ermsg( "SWACOS", 1, 0, "Deriv of ACOS(x) is infinite at x = -1 or +1" , '.' ); } else { ermsg( "SWASIN", 1, 0, "Deriv of ASIN(x) is infinite at x = -1 or +1" , '.' ); } return; } ns = n16 - 1; if (ns == 0) goto L_66; S1[2] = -2.0e0*U16[1]; if (ns == 1) goto L_66; S1[3] = -2.0e0; if (ns == 2) goto L_66; for (i = 3; i <= ns; i++) { S1[i + 1] = 0.0e0; } /* Now S1() = 1.0 - U**2 and derivs w.r.t. U. * Compute S2() = sqrt(1.0 - U**2) and derivs w.r.t. U. * Then change sign if doing Arc Cosine. * */ L_66: ; swsqrt( ns, s1, s2 ); if (acosin) swpro1( ns, -1.0e0, s2, s2 ); /* Compute 1.0/S2 and derivs w.r.t. U. Store result starting in Y(2) * */ swpwri( ns, -1, s2, &Y16[2] ); /* Convert Y from derivs w.r.t. U to derivs w.r.t. T */ swchn( n16, u16, y16 ); return; } /* end of function */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ void /*FUNCTION*/ swatn2( long n9, float y9[], float x9[], float a9[]) { long int i, nm1, np1; float big, s1[NMAXP1], s2[NMAXP1], s3[NMAXP1], s4[NMAXP1], xx, yy; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const A9 = &a9[0] - 1; float *const S1 = &s1[0] - 1; float *const S2 = &s2[0] - 1; float *const S3 = &s3[0] - 1; float *const S4 = &s4[0] - 1; float *const X9 = &x9[0] - 1; float *const Y9 = &y9[0] - 1; /* end of OFFSET VECTORS */ /* COMPUTES.. A = ATAN2(Y,X) = ATAN(Y/X) and derivs w.r.t. T * */ A9[1] = atan2f( Y9[1], X9[1] ); switch (IARITHIF(n9 - 1)) { case -1: goto L_90; case 0: goto L_92; case 1: goto L_94; } L_90: return; /* SPECIAL FOR N9 = 1 */ L_92: ; big = fmaxf( fabsf( X9[1] ), fabsf( Y9[1] ) ); yy = Y9[1]/big; xx = X9[1]/big; A9[2] = (xx*(Y9[2]/big) - yy*(X9[2]/big))/(SQ(xx) + SQ(yy)); return; /* GENERAL CASE N9 .ge. 2 */ L_94: ; np1 = n9 + 1; big = fmaxf( fabsf( X9[1] ), fabsf( Y9[1] ) ); for (i = 1; i <= np1; i++) { S1[i] = X9[i]/big; S2[i] = Y9[i]/big; } nm1 = n9 - 1; swpro( nm1, s1, &S2[2], s3 ); swpro( nm1, s2, &S1[2], s4 ); swdif( nm1, s3, s4, s3 ); swpro( nm1, s1, s1, s1 ); swpro( nm1, s2, s2, s2 ); swsum( nm1, s1, s2, s1 ); swquo( nm1, s3, s1, &A9[2] ); return; } /* end of function */ void /*FUNCTION*/ swsum( long nderiv, float u[], float v[], float w[]) { long int i, np1; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U = &u[0] - 1; float *const V = &v[0] - 1; float *const W = &w[0] - 1; /* end of OFFSET VECTORS */ /* COMPUTE.. W=U+V and derivs w.r.t. T * */ np1 = nderiv + 1; for (i = 1; i <= np1; i++) { W[i] = U[i] + V[i]; } return; } /* end of function */ void /*FUNCTION*/ swdif( long n2, float u2[], float v2[], float w2[]) { long int i, np1; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U2 = &u2[0] - 1; float *const V2 = &v2[0] - 1; float *const W2 = &w2[0] - 1; /* end of OFFSET VECTORS */ /* COMPUTE.. W2=U2-V2 and derivs w.r.t. T */ np1 = n2 + 1; for (i = 1; i <= np1; i++) { W2[i] = U2[i] - V2[i]; } return; } /* end of function */ void /*FUNCTION*/ swsqrt( long n15, float u15[], float y15[]) { long int i; float fac, s1[NMAXP1], scale, sqrtu; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const S1 = &s1[0] - 1; float *const U15 = &u15[0] - 1; float *const Y15 = &y15[0] - 1; /* end of OFFSET VECTORS */ /* COMPUTE.. Y= SQRT(U) and derivs w.r.t. T * */ sqrtu = sqrtf( U15[1] ); if (n15 == 0) goto L_152; if (sqrtu == 0.0e0) { ermsg( "SWSQRT", 2, 0, "Deriv of sqrtf(x) is infinite at x = 0" , '.' ); for (i = 2; i <= (n15 + 1); i++) { Y15[i] = 0.0e0; } return; } if (n15 > 1) goto L_154; Y15[2] = U15[2]/(sqrtu + sqrtu); L_152: Y15[1] = sqrtu; return; /* GENERAL CASE FOR N .ge. 2 */ L_154: ; scale = 1.0e0/U15[1]; Y15[1] = 1.0e0; S1[1] = 1.0e0; fac = 0.5e0; for (i = 1; i <= n15; i++) { Y15[i + 1] = Y15[i]*fac; fac -= 1.0e0; S1[i + 1] = U15[i + 1]*scale; } swchn( n15, s1, y15 ); Y15[1] = sqrtu; for (i = 1; i <= n15; i++) { Y15[i + 1] *= sqrtu; } return; } /* end of function */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * * */ void /*FUNCTION*/ swexp( long n4, float u4[], float y4[]) { long int i; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U4 = &u4[0] - 1; float *const Y4 = &y4[0] - 1; /* end of OFFSET VECTORS */ /* COMPUTE.. Y4=EXP(U4) and derivs w.r.t. T * */ Y4[1] = expf( U4[1] ); if (n4 == 0) return; for (i = 1; i <= n4; i++) { Y4[i + 1] = Y4[i]; } swchn( n4, u4, y4 ); return; } /* end of function */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ void /*FUNCTION*/ swsin( long n7, float u7[], float y7[]) { long int i; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U7 = &u7[0] - 1; float *const Y7 = &y7[0] - 1; /* end of OFFSET VECTORS */ /* COMPUTE.. Y7=SIN(U7) and derivs w.r.t. T * */ Y7[1] = sinf( U7[1] ); if (n7 == 0) return; Y7[2] = cosf( U7[1] ); if (n7 == 1) { Y7[2] *= U7[2]; } else { for (i = 2; i <= n7; i++) { Y7[i + 1] = -Y7[i - 1]; } swchn( n7, u7, y7 ); } return; } /* end of function */ void /*FUNCTION*/ swcos( long n7, float u7[], float y7[]) { long int i; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U7 = &u7[0] - 1; float *const Y7 = &y7[0] - 1; /* end of OFFSET VECTORS */ /* COMPUTE.. Y7=COS(U7) and derivs w.r.t. T * */ Y7[1] = cosf( U7[1] ); if (n7 == 0) return; Y7[2] = -sinf( U7[1] ); if (n7 == 1) { Y7[2] *= U7[2]; } else { for (i = 2; i <= n7; i++) { Y7[i + 1] = -Y7[i - 1]; } swchn( n7, u7, y7 ); } return; } /* end of function */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ void /*FUNCTION*/ swtan( long n7, float u7[], float y7[]) { long int i, n7m1; float cs, s1[NMAXP1], sn; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const S1 = &s1[0] - 1; float *const U7 = &u7[0] - 1; float *const Y7 = &y7[0] - 1; /* end of OFFSET VECTORS */ /* COMPUTE.. Y7=TAN(U7) and derivs w.r.t. T * The first deriv is 1/cos(u)**2 * */ if (n7 == 0) { Y7[1] = tanf( U7[1] ); return; } sn = sinf( U7[1] ); cs = cosf( U7[1] ); Y7[1] = sn/cs; /* Compute Cos(U) and its derivs w.r.t. U in S1(2..N+1) * */ S1[2] = cs; S1[3] = -sn; for (i = 4; i <= (n7 + 1); i++) { S1[i] = -S1[i - 2]; } n7m1 = n7 - 1; /* Convert to derivs w.r.t. t */ swchn( n7m1, u7, &S1[2] ); /* Compute Cos(U)**2 & derivs w.r.t. t */ swpro( n7m1, &S1[2], &S1[2], &S1[2] ); /* Divide first deriv of U by Cos(U)**2, * both with derivs w.r.t. t */ swquo( n7m1, &U7[2], &S1[2], &Y7[2] ); return; } /* end of function */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ void /*FUNCTION*/ swsinh( long n7, float u7[], float y7[]) { long int i; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U7 = &u7[0] - 1; float *const Y7 = &y7[0] - 1; /* end of OFFSET VECTORS */ /* COMPUTE.. Y7=SINH(U7) and derivs w.r.t. T * */ Y7[1] = sinhf( U7[1] ); if (n7 == 0) return; Y7[2] = coshf( U7[1] ); if (n7 == 1) { Y7[2] *= U7[2]; } else { for (i = 2; i <= n7; i++) { Y7[i + 1] = Y7[i - 1]; } swchn( n7, u7, y7 ); } return; } /* end of function */ void /*FUNCTION*/ swcosh( long n7, float u7[], float y7[]) { long int i; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U7 = &u7[0] - 1; float *const Y7 = &y7[0] - 1; /* end of OFFSET VECTORS */ /* COMPUTE.. Y7=COSH(U7) and derivs w.r.t. T * */ Y7[1] = coshf( U7[1] ); if (n7 == 0) return; Y7[2] = sinhf( U7[1] ); if (n7 == 1) { Y7[2] *= U7[2]; } else { for (i = 2; i <= n7; i++) { Y7[i + 1] = Y7[i - 1]; } swchn( n7, u7, y7 ); } return; } /* end of function */ void /*FUNCTION*/ swtanh( long n7, float u7[], float y7[]) { long int i, n7m1; float ch, s1[NMAXP1], sh; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const S1 = &s1[0] - 1; float *const U7 = &u7[0] - 1; float *const Y7 = &y7[0] - 1; /* end of OFFSET VECTORS */ /* COMPUTE.. Y7=TANH(U7) and derivs w.r.t. T * The first deriv is 1/cosh(u)**2 * */ if (n7 == 0) { Y7[1] = tanhf( U7[1] ); return; } sh = sinhf( U7[1] ); ch = coshf( U7[1] ); Y7[1] = sh/ch; /* Compute Cosh(U) and its derivs w.r.t. U in S1(2..N+1) * */ S1[2] = ch; S1[3] = sh; for (i = 4; i <= (n7 + 1); i++) { S1[i] = S1[i - 2]; } n7m1 = n7 - 1; /* Convert to derivs w.r.t. t */ swchn( n7m1, u7, &S1[2] ); /* Compute Cosh(U)**2 & derivs w.r.t. t */ swpro( n7m1, &S1[2], &S1[2], &S1[2] ); /* Divide first deriv of U by Cosh(U)**2, * both with derivs w.r.t. t */ swquo( n7m1, &U7[2], &S1[2], &Y7[2] ); return; } /* end of function */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ void /*FUNCTION*/ swset( long n10, float uval, float uder, float u10[]) { long int i; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U10 = &u10[0] - 1; /* end of OFFSET VECTORS */ /* INITIALIZE THE ARRAY U10().. SET VALUE=UVAL, 1-ST DERIV = UDER, * AND HIGHER derivs = ZERO */ switch (IARITHIF(n10 - 1)) { case -1: goto L_104; case 0: goto L_102; case 1: goto L_100; } L_100: for (i = 2; i <= n10; i++) { U10[i + 1] = 0.0e0; } L_102: U10[2] = uder; L_104: U10[1] = uval; return; } /* end of function */ void /*FUNCTION*/ swsum1( long n11, float c11, float u11[], float y11[]) { long int i; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U11 = &u11[0] - 1; float *const Y11 = &y11[0] - 1; /* end of OFFSET VECTORS */ /* COMPUTE.. Y=C+U WHERE C IS A SCALAR * */ Y11[1] = c11 + U11[1]; for (i = 1; i <= n11; i++) { Y11[i + 1] = U11[i + 1]; } return; } /* end of function */ void /*FUNCTION*/ swdif1( long n12, float c12, float u12[], float y12[]) { long int i; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U12 = &u12[0] - 1; float *const Y12 = &y12[0] - 1; /* end of OFFSET VECTORS */ /* COMPUTE.. Y=C-U WHERE C IS A SCALAR * */ Y12[1] = c12 - U12[1]; for (i = 1; i <= n12; i++) { Y12[i + 1] = -U12[i + 1]; } return; } /* end of function */ void /*FUNCTION*/ swpro1( long n13, float c13, float u13[], float y13[]) { long int i; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U13 = &u13[0] - 1; float *const Y13 = &y13[0] - 1; /* end of OFFSET VECTORS */ /* COMPUTE.. Y=C*U WHERE C IS A SCALAR * */ for (i = 1; i <= (n13 + 1); i++) { Y13[i] = c13*U13[i]; } return; } /* end of function */ void /*FUNCTION*/ swquo1( long n14, float c14, float u14[], float y14[]) { long int i; float s1[NMAXP1]; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const S1 = &s1[0] - 1; float *const U14 = &u14[0] - 1; float *const Y14 = &y14[0] - 1; /* end of OFFSET VECTORS */ /* COMPUTE.. Y=C/U WHERE C IS A SCALAR * */ if (U14[1] == 0.0e0) { ermsg( "SWQUO1", 3, 0, "Zero divisor.", '.' ); return; } S1[1] = c14; for (i = 2; i <= (n14 + 1); i++) { S1[i] = 0.0e0; } swquo( n14, s1, u14, y14 ); return; } /* end of function */ void /*FUNCTION*/ swlog( long nderiv, float u[], float y[]) { long int ii; float fac, uinv; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U = &u[0] - 1; float *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* C. L. LAWSON, JPL, 1971 NOV 19 * * COMPUTE.. Y= LOG BASE E OF U and derivs w.r.t. T * ERROR STOP IN this subr IF U(1) .le. 0. * */ Y[1] = logf( U[1] ); if (nderiv <= 1) { if (nderiv == 1) Y[2] = U[2]/U[1]; return; } uinv = 1.0e0/U[1]; Y[2] = uinv; fac = -uinv; for (ii = 2; ii <= nderiv; ii++) { Y[ii + 1] = Y[ii]*fac; fac -= uinv; } /* Y() NOW CONTAINS VALUE and derivs w.r.t. U * USE CHAIN RULE TO convert TO derivs w.r.t. T * */ swchn( nderiv, u, y ); return; } /* end of function */ void /*FUNCTION*/ swpwri( long nderiv, long mpwr, float u[], float y[]) { LOGICAL32 short_; long int i, ii, ii1, ii2, j, m; float d, fac, uinv; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U = &u[0] - 1; float *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* COMPUTE.. Y = U**MPWR and derivs w.r.t. T * MPWR is an integer, POS.,NEG.,OR ZERO. * MPWR does not depend ON T. * U may depend on T * Issues an error message if both (U(1) .eq. 0.) and (MPWR .lt. 0) * If MPWR .eq. 0 then set Y(1) = 1. and all derivs = 0. * */ ii2 = nderiv; short_ = FALSE; m = mpwr; switch (IARITHIF(m)) { case -1: goto L_110; case 0: goto L_35; case 1: goto L_45; } /* HERE MPWR .eq. 0 */ L_35: Y[1] = 1.0e0; if (nderiv <= 0) return; for (j = 1; j <= nderiv; j++) { Y[j + 1] = 0.0e0; } return; /* HERE MPWR IS POSITIVE */ L_45: if (U[1] != 0.0e0) goto L_50; Y[1] = 0.0e0; if (nderiv == 0) return; if (m <= nderiv) goto L_46; ii = 1; goto L_130; L_46: for (i = 1; i <= nderiv; i++) { Y[i + 1] = 0.0e0; } fac = 1.0e0; d = 0.0e0; for (i = 1; i <= m; i++) { d += 1.0e0; fac *= d; } Y[m + 1] = fac; goto L_150; L_50: if (m >= nderiv) goto L_112; ii2 = m; short_ = TRUE; goto L_112; /* HERE MPWR .lt 0 */ L_110: ; if (U[1] == 0.0e0) { ierm1( "SWPWRI", 4, 0, "U**M is infinite when U = 0. and M < 0" , "M", m, '.' ); for (i = 1; i <= (nderiv + 1); i++) { Y[i] = 0.0e0; } return; } L_112: ; Y[1] = powif(U[1],m); if (nderiv == 0) return; uinv = 1.0e0/U[1]; fac = uinv*m; ii1 = 1; for (ii = ii1; ii <= ii2; ii++) { Y[ii + 1] = Y[ii]*fac; fac -= uinv; } if (!short_) goto L_150; ii = ii2 + 1; /* SET HIGHER derivs TO ZERO. * */ L_130: ; for (j = ii; j <= nderiv; j++) { Y[j + 1] = 0.0e0; } /* Y() NOW CONTAINS VALUE and derivs w.r.t. U * USE CHAIN RULE TO convert TO derivs w.r.t. T * */ L_150: swchn( nderiv, u, y ); return; } /* end of function */ void /*FUNCTION*/ swchn( long nderiv, float u[], float y[]) { long int i, l, np1; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U = &u[0] - 1; float *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* Implements the chain rule. * Given values of y(u) and its derivs with respect to u in Y(), and * u(x) and its derivs with respect to x in U(), this subr replaces * the contents of Y() with y(u(x)) and its derivs with respect to x. * * NDERIV HIGHEST ORDER DERIVATIVE * * (U(I),I=1,NDERIV+1) INPUT.. value of u and derivs w.r.t. x * * (Y(I),I=1,NDERIV+1) INPUT.. value of y and derivs w.r.t. u * OUTPUT.. (Y(I+1),I=1,NDERIV) will be * replaced by derivs w.r.t. x */ if (nderiv == 0) return; if (U[2] != 1.) goto L_20; /* U(2) .eq. 1., See if higher derivs are 0. */ if (nderiv == 1) return; for (i = 2; i <= nderiv; i++) { if (U[i + 1] != 0.) goto L_50; } return; /* Test if all derivs of U are 0. */ L_20: for (i = 1; i <= nderiv; i++) { if (U[i + 1] != 0.) goto L_50; } for (i = 1; i <= nderiv; i++) { Y[i + 1] = 0.; } return; /* Here for the general case when U is not T and not constant. * */ L_50: np1 = nderiv + 1; for (l = 0; l <= (nderiv - 1); l++) { swpro( l, &Y[np1 - l], &U[2], &Y[np1 - l] ); } return; } /* end of function */ void /*FUNCTION*/ swrchn( long nderiv, float u[], float y[]) { long int l, np1; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U = &u[0] - 1; float *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* Implements the "reverse" chain rule. * Given values of y(u(x)) and its derivs with respect to x in Y(), * and u(x) and its derivs with respect to x in U(), this subr replac * the contents of Y() with y(u) and its derivs with respect to u. * Requires u'(x) .ne. 0. * * Note that this subr can be used to compute a representation of * the inverse function of a given function. This transformation * is called "reversion" of a power series. To do this let y() be * the inverse function of u() in a neighborhood of a point, x0. * Thus y(u(x)) = x for all x in a neighborhood of x0. * Then the value of y() and its derivs with respect to x, evaluated * at x0 are (x0, 1.0, 0.0, ..., 0.0). * If Y() contains these values on entry, and U() contains u(x) and * its derivs w.r.t. x, evaluated at x0, then on return Y() will * contain y(u) and its derivs w.r.t. u, evaluated at u0 = u(x0). * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * NDERIV [in] Highest order derivative to be considered. * * (U(I),I=1,NDERIV+1) [in] Contains the value of u(x) and its * first NDERIV derivatives with respect to x, evaluated at x0. * Require U(2) .ne. 0. * * (Y(I),I=1,NDERIV+1) [inout] On entry, Y() must contain y(u(x)) * and its first NDERIV derivatives w.r.t. x, evaluated at x0. * On return, Y() will contain y(u) and its first NDERIV * derivatives w.r.t. u, evaluated at u0 = u(x0). * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ np1 = nderiv + 1; if (U[2] == 0.0e0 && np1 >= 2) { ermsg( "SWRCHN", 9, 0, "U(2) is zero.", '.' ); return; } for (l = 2; l <= np1; l++) { swquo( np1 - l, &Y[l], &U[2], &Y[l] ); } return; } /* end of function */ /* PARAMETER translations */ #define NCOEF (1 + (NMAX*(NMAX + 1))/2) /* end of PARAMETER translations */ /* COMMON translations */ struct t_swcom { float c[NCOEF]; } swcom; /* end of COMMON translations */ void /*FUNCTION*/ swpro( long nderiv, float u[], float v[], float w[]) { long int i, id2, imid, ip1, j, jc, np1; float temp; static LOGICAL32 first = TRUE; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const C = &swcom.c[0] - 1; float *const U = &u[0] - 1; float *const V = &v[0] - 1; float *const W = &w[0] - 1; /* end of OFFSET VECTORS */ /* COMPUTE W(1)=U(1)*V(1) AND DERIVATIVES . * NDERIV HIGHEST ORDER DERIVATIVE TO BE COMPUTED * NDERIV .ge. 0 * (U(I),I=1,NDERIV+1) INPUT ARRAY.. * U(1)= VALUE * U(I+1)= I-TH DERIVATIVE OF U(1) * (V(I),I=1,NDERIV+1) INPUT ARRAY.. SAME FORMAT AS U(). * * (W(I),I=1,NDERIV+1) OUTPUT ARRAY.. * W(1)= U(1)*V(1) * W(I+1)=Ith deriv of W(1) * It is permissible for W() to occupy the same storage locations as * U() and/or V(). * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * NMAX Setting NMAX in the parameter statement determines * the highest order derivative this program unit will be able * to handle. The dimension of the internal saved array C() * will be set as a function of NMAX. * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ if (first) { spascl( NMAX, swcom.c ); first = FALSE; } if (nderiv > NMAX) { ierm1( "SWPRO", 5, 0, "Dimension NMAX in SWATAN, SWSUM, & SWPRO must be .ge. NDERIV" , "NMAX", NMAX, ',' ); ierv1( "NDERIV", nderiv, '.' ); return; } np1 = nderiv + 1; /* JC = (NDERIV*(NDERIV-1))/2 * DO 20 I=NP1,1,-1 * TEMP=0.0E0 * IP1=I+1 * DO 10 J=1,I * TEMP = TEMP + U(J)*V(IP1-J)*C(JC+J) * 10 continue * W(I) = TEMP * JC=JC- I+2 * 20 continue * * The following code does the same as the above commented-out * code but reduces the number of multiplies. * Treats the cases of NDERIV = 0, 1, 2, 3, and 4 special for * efficiency, since it is expected that most usage will involve * small values of NDERIV. In particular the chain rule subroutine, * _WCHN, will call this subr with a sequence of values of NDERIV * going down to 0. * Treats NDERIV > 4 in a general way, but takes account of * the fact that the C() terms in the "do 10" loop above are * symmetric about the middle term and the first and last term * are each = 1. * */ switch (np1) { case 1: goto L_41; case 2: goto L_42; case 3: goto L_43; case 4: goto L_44; case 5: goto L_45; } /* Fortran 77 specifies that a Computed Go To drops through is the * index is out of range. Thus we arrive here if NP1 < 1 or > 5. * NP1 < 1 would be an error. We choose not to use time to test for * it however. * NP1 > 5 is valid. In that case we do the following loop for I * running down from NP1 to 6, and then do the special statements * for 5, 4, 3, 2, and 1. * */ jc = (nderiv*(nderiv - 1))/2; for (i = np1; i >= 6; i--) { temp = U[1]*V[i] + U[i]*V[1]; ip1 = i + 1; id2 = i/2; for (j = 2; j <= id2; j++) { temp += C[jc + j]*(U[j]*V[ip1 - j] + U[ip1 - j]*V[j]); } if (2*id2 == i) { W[i] = temp; } else { imid = id2 + 1; W[i] = temp + C[jc + imid]*U[imid]*V[imid]; } jc += -i + 2; } L_45: W[5] = U[1]*V[5] + 4.0e0*(U[2]*V[4] + U[4]*V[2]) + 6.0e0*U[3]* V[3] + U[5]*V[1]; L_44: W[4] = U[1]*V[4] + 3.0e0*(U[2]*V[3] + U[3]*V[2]) + U[4]*V[1]; L_43: W[3] = U[1]*V[3] + 2.0e0*U[2]*V[2] + U[3]*V[1]; L_42: W[2] = U[1]*V[2] + U[2]*V[1]; L_41: W[1] = U[1]*V[1]; return; } /* end of function */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ void /*FUNCTION*/ swquo( long nderiv, float u[], float v[], float w[]) { long int i, id2, imid, ip1, j, jc, np1; float fac, temp; static LOGICAL32 first = TRUE; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const C = &swcom.c[0] - 1; float *const U = &u[0] - 1; float *const V = &v[0] - 1; float *const W = &w[0] - 1; /* end of OFFSET VECTORS */ /* COMPUTE.. W = U / V and derivs w.r.t. T * * V() must be distinct in storage from W(), but U() can be the same * array as W(). * * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ if (first) { spascl( NMAX, swcom.c ); first = FALSE; } if (nderiv > NMAX) { ierm1( "SWQUO", 6, 0, "Dimension NMAX in SWATAN, SWSUM, & SWPRO must be .ge. NDERIV" , "NMAX", NMAX, ',' ); ierv1( "NDERIV", nderiv, '.' ); return; } if (V[1] == 0.0e0) { ermsg( "SWQUO", 7, 0, "Zero divisor.", '.' ); return; } fac = 1.0e0/V[1]; W[1] = U[1]*fac; np1 = nderiv + 1; jc = 0; /* DO 40 I=2,NP1 * TEMP = U(I) * IP1=I+1 * DO 30 J=2,I * TEMP = TEMP - V(J)*W(IP1-J)*C(JC+J) * 30 CONTINUE * W(I) = TEMP * FAC * JC=JC + I-1 * 40 CONTINUE */ for (i = 2; i <= np1; i++) { temp = U[i] - W[1]*V[i]; ip1 = i + 1; id2 = i/2; for (j = 2; j <= id2; j++) { temp += -C[jc + j]*(W[j]*V[ip1 - j] + W[ip1 - j]*V[j]); } if (2*id2 != i) { imid = id2 + 1; temp += -C[jc + imid]*W[imid]*V[imid]; } W[i] = temp*fac; jc += i - 1; } return; } /* end of function */ void /*FUNCTION*/ spascl( long n, float c[]) { long int i, j, k, l; static LOGICAL32 done = FALSE; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const C = &c[0] - 1; /* end of OFFSET VECTORS */ /* C.L.LAWSON,JPL, 1969 DEC 10 * * N MAXIMUM ORDER DERIVATIVE TO BE * COMPUTED USING W-ARITHMETIC * * (C(I),I=1,NC) FIRST N+1 ROWS OF PASCAL TRIANGLE PACKED * WITH 1-ST DIAGONAL OMITTED AFTER 1-ST ROW. * NC= (N*(N+1)/2) + 1 * EXAMPLE.. IF N=3 then TRIANGLE IS 1 * 1 1 * 1 2 1 * 1 3 3 1 * * WHICH WILL BE STORED AS.. 1,1,2,1,3,3,1 * */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ if (done) return; done = TRUE; if (n < 0) { ierm1( "SPASCL", 8, 0, "Require N .ge. 0", "N", n, '.' ); return; } C[1] = 1.0e0; if (n == 0) return; C[2] = 1.0e0; if (n == 1) return; k = 3; for (i = 1; i <= (n - 1); i++) { l = i + 1; for (j = 1; j <= i; j++) { C[k] = C[k - l] + C[k - l + 1]; k += 1; } C[k] = 1.0e0; k += 1; } return; } /* end of function */