/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:44 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "ducomp.h" #include /* COMMON translations */ struct t_ucom1 { long int n, m1, m2; } ucom1; struct t_ucom2 { long int l1, l2; } ucom2; /* end of COMMON translations */ void /*FUNCTION*/ dusetn( long nin, long m1in, long m2in) { /* Base name of file: DUCOMP * Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 2001-06-18 DUCOMP Krogh Replaced "M1 .eq.. 0" with "M1 .eq. 0" *>> 1996-04-30 DUCOMP Krogh Removed redundant save statement. *>> 1994-11-02 DUCOMP Krogh Changes to use M77CON *>> 1994-08-04 DUCOMP CLL New subrs: [D|S]USETN & [D|S]UGETN. *>> 1993-11-11 CLL Changed Entries to separate Subroutines. *>> 1990-01-23 CLL Replace M with MPWR in call to IERM1 *>> 1987-12-07 DUCOMP Lawson Initial code. * * The file DUCOMP contains a set of program units to perform * computation on arrays representing the value of a quantity * and (optionally) its first and second partial derivatives * with respect to a set of N independent variables. * * The lowest and highest order derivs desired are specified by * the integers M1 and M2 which must satisfy * 0 .le. M1 .le. M2 .le. 2 * The integers N, M1, and M2 must be set by the user by a call * to [D|S]USETN, which will put these values in the COMMON block UCOM1. * * Each array containing a U-variable for which second partial * derivs are to be carried must be dimensioned at least * ((N+2)*(N+1))/2 EXAMPLES.. * N = 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 * DIMENSION = 3 6 10 15 21 28 36 45 55 66 78 91 105 120 136 153 * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Subroutines * * subroutine DUSETN ( N, M1, M2) * subroutine DUGETN ( N, M1, M2, L1, L2) * subroutine DUSET ( VAL, KEY, Y9) * * subroutine DUPRO ( U,V,Y) * subroutine DUQUO ( U1,V1,Y1) * subroutine DUSUM ( U3,V3,Y3) * subroutine DUDIF ( U4,V4,Y4) * subroutine DUSUM1 ( C5,V5,Y5) * subroutine DUDIF1 ( C6,V6,Y6) * subroutine DUPRO1 ( C7,V7,Y7) * subroutine DUQUO1 ( C8,V8,Y8) * * subroutine DUSQRT (U,Y) * subroutine DUEXP (U1,Y1) * subroutine DULOG (U2,Y2) * subroutine DUPWRI (MPWR,U3,Y3) * * subroutine DUSIN (U,Y) * subroutine DUCOS (U ,Y ) * subroutine DUSINH (U ,Y ) * subroutine DUCOSH (U ,Y ) * subroutine DUATAN (U1,Y1) * subroutine DUATN2 (V2,U2,Y2) * subroutine DUASIN (U, Y) * subroutine DUACOS (U, Y) * subroutine DUTAN (U, Y) * subroutine DUTANH (U, Y) * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *--D replaces "?": ?UCOMP, ?UACOS, ?UASIN, ?UATAN, ?UATN2, ?UCOS *-- & ?UCOSH, ?UDIF, ?UDIF1, ?UEXP, ?UGETN, ?ULOG, ?UPRO, ?UPRO1 *-- & ?UPWRI, ?UQUO, ?UQUO1, ?USET, ?USETN, ?USIN, ?USINH, ?USQRT *-- & ?USUM, ?USUM1, ?UTAN, ?UTANH, ?UACS * Subprograms referecnced by both versions: ERMSG * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * C. L. Lawson, JPL, 1969 Dec 4 * 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, Sept 1987. * ================================================================== * subroutine DUSETN ( Nin, M1in, M2in) * * Nin [in] Specifies the number of independent variables. * Require Nin .ge. 1. * M1in, M2in [in] These specify the lowest and highest order derivs * desired. These must satisfy 0 .le. M1in .le. M2in .le. 2. * * This subr copies Nin, M1in, and M2in into common variables N, M1, and * M2. * It also computes L1 and L2 based on N, M1, and M2, and stores * L1 and L2 into common. L1 and L2 specify the first and last * locations in a U-variable array that will be operated on when * dealing with derivs of orders M1 through M2 and * N independent variables. * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ucom1.n = nin; ucom1.m1 = m1in; ucom1.m2 = m2in; /* SET L1 AND L2 */ if (ucom1.m1 == 0) { ucom2.l1 = 1; } else if (ucom1.m1 == 1) { ucom2.l1 = 2; } else { ucom2.l1 = ucom1.n + 2; } if (ucom1.m2 == 0) { ucom2.l2 = 1; } else if (ucom1.m2 == 1) { ucom2.l2 = ucom1.n + 1; } else { ucom2.l2 = 1 + ucom1.n + ((ucom1.n*(ucom1.n + 1))/2); } return; } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dugetn( long *nout, long *m1out, long *m2out, long *l1out, long *l2out) { /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ *nout = ucom1.n; *m1out = ucom1.m1; *m2out = ucom1.m2; *l1out = ucom2.l1; *l2out = ucom2.l2; return; } /* end of function */ /* ================================================================== */ /* PARAMETER translations */ #define ONE 1.0e0 #define ZERO 0.0e0 /* end of PARAMETER translations */ void /*FUNCTION*/ duset( double val, long key, double y9[]) { long int i, j, k; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Y9 = &y9[0] - 1; /* end of OFFSET VECTORS */ /* A convenience routine for assigning an initial value to the * U-variable, Y9(). The parts of Y9() assigned will * depend on M1 and M2. * * If M1 .le. 0 .le. M2, Y9() will be assigned the value, VAL. * * If M1 .le. 1 .le. M2, the first partial derivs of Y9() will be * assigned depending on KEY. KEY must satisfy 0 .le. KEY .le. N. * If KEY .eq.. 0, the U-variable is to be constant relative to the N * independent variables. Thus all first partial derivs will be set * to zero. * If KEY is in [1,N], Y9() is being set as th KEYth independent * variable. Thus the KEYth first partial deriv will be set to 1.0 * and all other first partial derivs will be set to zero. * * If M1 .le. 2 .le. M2, the second partial derivs of Y9() will be * set to zero. * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ switch (IARITHIF(ucom1.m1 - 1)) { case -1: goto L_230; case 0: goto L_240; case 1: goto L_250; } /* VALUE */ L_230: Y9[1] = val; if (ucom1.m2 == 0) return; /* 1ST PARTIALS */ L_240: ; for (i = 1; i <= ucom1.n; i++) { Y9[i + 1] = ZERO; } if (key >= 1 && key <= ucom1.n) Y9[key + 1] = ONE; if (ucom1.m2 == 1) return; /* 2ND PARTIALS */ L_250: ; k = ucom1.n + 2; for (i = 1; i <= ucom1.n; i++) { for (j = 1; j <= i; j++) { Y9[k] = ZERO; k += 1; } } return; } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dupro( double u[], double v[], double y[]) { long int i, j, k; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const U = &u[0] - 1; double *const V = &v[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* Computes Y = U * V * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ switch (IARITHIF(ucom1.m2 - 1)) { case -1: goto L_60; case 0: goto L_40; case 1: goto L_10; } /* 2ND PARTIALS */ L_10: ; k = ucom1.n + 2; for (i = 1; i <= ucom1.n; i++) { for (j = 1; j <= i; j++) { Y[k] = U[1]*V[k] + U[k]*V[1] + U[i + 1]*V[j + 1] + U[j + 1]* V[i + 1]; k += 1; } } if (ucom1.m1 == 2) return; /* 1ST PARTIALS */ L_40: ; for (i = 1; i <= ucom1.n; i++) { Y[i + 1] = U[1]*V[i + 1] + U[i + 1]*V[1]; } if (ucom1.m1 == 1) return; /* VALUE */ L_60: ; Y[1] = U[1]*V[1]; return; } /* end of function */ void /*FUNCTION*/ duquo( double u1[], double v1[], double y1[]) { long int i, j, k; double vinv; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const U1 = &u1[0] - 1; double *const V1 = &v1[0] - 1; double *const Y1 = &y1[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* COMPUTE.. Y=U/V * */ if (V1[1] == ZERO) { ermsg( "DUQUO", 5, 0, "Denominator is zero.", '.' ); return; } switch (IARITHIF(ucom1.m1 - 1)) { case -1: goto L_80; case 0: goto L_90; case 1: goto L_100; } /* VALUE */ L_80: ; Y1[1] = U1[1]/V1[1]; if (ucom1.m2 == 0) return; /* 1ST PARTIAL DERIVS */ L_90: ; vinv = ONE/V1[1]; for (i = 1; i <= ucom1.n; i++) { Y1[i + 1] = (U1[i + 1] - Y1[1]*V1[i + 1])*vinv; } if (ucom1.m2 == 1) return; goto L_110; /* 2ND PARTIAL DERIVS */ L_100: vinv = ONE/V1[1]; L_110: ; k = ucom1.n + 2; for (i = 1; i <= ucom1.n; i++) { for (j = 1; j <= i; j++) { Y1[k] = (U1[k] - Y1[1]*V1[k] - Y1[i + 1]*V1[j + 1] - Y1[j + 1]* V1[i + 1])*vinv; k += 1; } } return; } /* end of function */ void /*FUNCTION*/ dusum( double u3[], double v3[], double y3[]) { long int k; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const U3 = &u3[0] - 1; double *const V3 = &v3[0] - 1; double *const Y3 = &y3[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* COMPUTE Y=U+V AND PARTIAL DERIVS * */ for (k = ucom2.l1; k <= ucom2.l2; k++) { Y3[k] = U3[k] + V3[k]; } return; } /* end of function */ void /*FUNCTION*/ dudif( double u4[], double v4[], double y4[]) { long int k; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const U4 = &u4[0] - 1; double *const V4 = &v4[0] - 1; double *const Y4 = &y4[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* COMPUTE Y=U-V AND PARTIAL DERIVS * */ for (k = ucom2.l1; k <= ucom2.l2; k++) { Y4[k] = U4[k] - V4[k]; } return; } /* end of function */ void /*FUNCTION*/ dusum1( double c5, double v5[], double y5[]) { long int k; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const V5 = &v5[0] - 1; double *const Y5 = &y5[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* COMPUTE Y=C+V WHERE C IS A CONSTANT * */ for (k = ucom2.l1; k <= ucom2.l2; k++) { Y5[k] = V5[k]; } if (ucom1.m1 == 0) Y5[1] += c5; return; } /* end of function */ void /*FUNCTION*/ dudif1( double c6, double v6[], double y6[]) { long int k; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const V6 = &v6[0] - 1; double *const Y6 = &y6[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* COMPUTE Y=C-V WHERE C IS A CONSTANT * */ for (k = ucom2.l1; k <= ucom2.l2; k++) { Y6[k] = -V6[k]; } /* Using + in following statement because have just set * Y6(1) = -V6(1). * */ if (ucom1.m1 == 0) Y6[1] += c6; return; } /* end of function */ void /*FUNCTION*/ dupro1( double c7, double v7[], double y7[]) { long int k; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const V7 = &v7[0] - 1; double *const Y7 = &y7[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* COMPUTE Y=C*V WHERE C IS A CONSTANT * */ for (k = ucom2.l1; k <= ucom2.l2; k++) { Y7[k] = c7*V7[k]; } return; } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ duquo1( double c8, double v8[], double y8[]) { long int i, j, k; double fac; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const V8 = &v8[0] - 1; double *const Y8 = &y8[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* COMPUTE Y=C/V WHERE C IS A CONSTANT * */ if (V8[1] == ZERO) { ermsg( "DUQUO1", 1, 0, "Denominator is zero.", '.' ); return; } switch (IARITHIF(ucom1.m1 - 1)) { case -1: goto L_204; case 0: goto L_208; case 1: goto L_214; } /* VALUE */ L_204: Y8[1] = c8/V8[1]; if (ucom1.m2 == 0) return; /* 1ST PARTIALS */ L_208: ; fac = -Y8[1]/V8[1]; for (i = 1; i <= ucom1.n; i++) { Y8[i + 1] = V8[i + 1]*fac; } if (ucom1.m2 == 1) return; /* 2ND PARTIALS */ L_214: fac = -ONE/V8[1]; k = ucom1.n + 2; for (i = 1; i <= ucom1.n; i++) { for (j = 1; j <= i; j++) { Y8[k] = (Y8[1]*V8[k] + Y8[i + 1]*V8[j + 1] + Y8[j + 1]* V8[i + 1])*fac; k += 1; } } return; } /* end of function */ /* ------------------------------------------------------------------ */ /* PARAMETER translations */ #define HALF 0.5e0 /* end of PARAMETER translations */ void /*FUNCTION*/ dusqrt( double u[], double y[]) { long int i, j, k; double d1, fac, fac2; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const U = &u[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* COMPUTE.. Y=SQRT(U) * * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ switch (IARITHIF(ucom1.m1 - 1)) { case -1: goto L_4; case 0: goto L_6; case 1: goto L_16; } L_4: Y[1] = sqrt( U[1] ); if (ucom1.m2 == 0) return; L_6: if (Y[1] == ZERO) { ermsg( "DUSQRT", 2, 0, "Deriv of sqrt(x) is infinite at x = 0" , '.' ); return; } d1 = HALF/Y[1]; for (i = 1; i <= ucom1.n; i++) { Y[i + 1] = U[i + 1]*d1; } if (ucom1.m2 == 1) return; goto L_18; L_16: ; d1 = HALF/Y[1]; L_18: k = ucom1.n + 1; fac = ONE/Y[1]; for (i = 1; i <= ucom1.n; i++) { fac2 = fac*Y[i + 1]; for (j = 1; j <= i; j++) { k += 1; Y[k] = d1*U[k] - Y[j + 1]*fac2; } } return; } /* end of function */ void /*FUNCTION*/ duexp( double u1[], double y1[]) { long int i, j, k; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const U1 = &u1[0] - 1; double *const Y1 = &y1[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* Y=EXP(U) */ switch (IARITHIF(ucom1.m1 - 1)) { case -1: goto L_30; case 0: goto L_31; case 1: goto L_34; } L_30: Y1[1] = exp( U1[1] ); if (ucom1.m2 == 0) return; L_31: ; for (i = 1; i <= ucom1.n; i++) { Y1[i + 1] = Y1[1]*U1[i + 1]; } if (ucom1.m2 == 1) return; L_34: ; k = ucom1.n + 1; for (i = 1; i <= ucom1.n; i++) { for (j = 1; j <= i; j++) { k += 1; Y1[k] = Y1[1]*(U1[k] + U1[i + 1]*U1[j + 1]); } } return; } /* end of function */ void /*FUNCTION*/ dulog( double u2[], double y2[]) { long int i, j, k; double d1; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const U2 = &u2[0] - 1; double *const Y2 = &y2[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* Y=LOG(U) */ switch (IARITHIF(ucom1.m1 - 1)) { case -1: goto L_40; case 0: goto L_42; case 1: goto L_48; } L_40: Y2[1] = log( U2[1] ); if (ucom1.m2 == 0) return; L_42: d1 = ONE/U2[1]; for (i = 1; i <= ucom1.n; i++) { Y2[i + 1] = U2[i + 1]*d1; } if (ucom1.m2 == 1) return; goto L_49; L_48: d1 = ONE/U2[1]; L_49: ; k = ucom1.n + 1; for (i = 1; i <= ucom1.n; i++) { for (j = 1; j <= i; j++) { k += 1; Y2[k] = d1*U2[k] - Y2[i + 1]*Y2[j + 1]; } } return; } /* end of function */ /* PARAMETER translations */ #define TWO 2.0e0 /* end of PARAMETER translations */ void /*FUNCTION*/ dupwri( long mpwr, double u3[], double y3[]) { long int i, j, k, k1, k2; double fac, fac2, fm; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const U3 = &u3[0] - 1; double *const Y3 = &y3[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* Y=U**MPWR MPWR IS A CONSTANT * MPWR MAY BE POS.,ZERO, OR NEG. IF MPWR = 0 THEN Y=1. INDEPENDENT * OF U. IF MPWR IS NEG. THEN ERROR STOP IF U IS ZERO. * */ if (mpwr == 0) goto L_56; if (U3[1] != ZERO) goto L_100; /* U = 0. OR MPWR = 0 */ if (mpwr < 0) { ierm1( "DUPWRI", 4, 0, "U**M is infinite when U = 0. and M < 0" , "M", mpwr, '.' ); return; } /* U=0. AND MPWR .GE. 0 */ L_56: switch (IARITHIF(ucom1.m1 - 1)) { case -1: goto L_64; case 0: goto L_68; case 1: goto L_84; } L_64: Y3[1] = ZERO; if (mpwr == 0) Y3[1] = ONE; if (ucom1.m2 == 0) return; L_68: if (mpwr == 1) goto L_74; for (i = 1; i <= ucom1.n; i++) { Y3[i + 1] = ZERO; } goto L_80; L_74: for (i = 1; i <= ucom1.n; i++) { Y3[i + 1] = U3[i + 1]; } L_80: if (ucom1.m2 == 1) return; L_84: ; if (mpwr == 2) { k = ucom1.n + 2; for (i = 1; i <= ucom1.n; i++) { for (j = 1; j <= i; j++) { Y3[k] = TWO*U3[i + 1]*U3[j + 1]; k += 1; } } return; } k1 = ucom1.n + 1; k2 = k1 - 1 + (ucom1.n*(ucom1.n + 1))/2; if (mpwr == 1) { for (k = k1; k <= k2; k++) { Y3[k] = U3[k]; } } else { for (k = k1; k <= k2; k++) { Y3[k] = ZERO; } } return; /* END.. U .eq. 0. .OR. MPWR .eq. 0 * * BEGIN.. U .NE. 0. .AND. MPWR .NE.0 */ L_100: ; switch (IARITHIF(ucom1.m1 - 1)) { case -1: goto L_104; case 0: goto L_108; case 1: goto L_112; } L_104: Y3[1] = powi(U3[1],mpwr); if (ucom1.m2 == 0) return; L_108: fm = mpwr; fac = fm*(Y3[1]/U3[1]); for (i = 1; i <= ucom1.n; i++) { Y3[i + 1] = fac*U3[i + 1]; } if (ucom1.m2 == 1) return; goto L_114; L_112: fm = mpwr; fac = fm*(Y3[1]/U3[1]); L_114: fac2 = (fm - ONE)/U3[1]; k = ucom1.n + 2; for (i = 1; i <= ucom1.n; i++) { for (j = 1; j <= i; j++) { Y3[k] = fac*(fac2*U3[i + 1]*U3[j + 1] + U3[k]); k += 1; } } return; /* END.. U .NE. 0. .AND. MPWR .NE. 0 */ } /* end of function */ /* ------------------------------------------------------------------ */ void /*FUNCTION*/ dusin( double u[], double y[]) { long int i, j, k; double d1, d2, fac; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const U = &u[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* computes Y=SIN(U) * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ if (ucom1.m1 <= 0) { Y[1] = sin( U[1] ); if (ucom1.m2 == 0) return; } d2 = -Y[1]; d1 = cos( U[1] ); if (ucom1.m1 <= 1) goto L_16; goto L_24; /* 1ST PARTIALS */ L_16: for (i = 1; i <= ucom1.n; i++) { Y[i + 1] = d1*U[i + 1]; } if (ucom1.m2 == 1) return; /* 2ND PARTIALS */ L_24: k = ucom1.n + 2; for (i = 1; i <= ucom1.n; i++) { fac = d2*U[i + 1]; for (j = 1; j <= i; j++) { Y[k] = fac*U[j + 1] + d1*U[k]; k += 1; } } return; } /* end of function */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ void /*FUNCTION*/ ducos( double u[], double y[]) { long int i, j, k; double d1, d2, fac; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const U = &u[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* COMPUTE.. Y=COS(U) * */ if (ucom1.m1 <= 0) { Y[1] = cos( U[1] ); if (ucom1.m2 == 0) return; } d2 = -Y[1]; d1 = -sin( U[1] ); if (ucom1.m1 <= 1) goto L_16; goto L_24; /* 1ST PARTIALS */ L_16: for (i = 1; i <= ucom1.n; i++) { Y[i + 1] = d1*U[i + 1]; } if (ucom1.m2 == 1) return; /* 2ND PARTIALS */ L_24: k = ucom1.n + 2; for (i = 1; i <= ucom1.n; i++) { fac = d2*U[i + 1]; for (j = 1; j <= i; j++) { Y[k] = fac*U[j + 1] + d1*U[k]; k += 1; } } return; } /* end of function */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ void /*FUNCTION*/ dusinh( double u[], double y[]) { long int i, j, k; double d1, d2, fac; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const U = &u[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* COMPUTE.. Y=SINH(U) * */ if (ucom1.m1 <= 0) { Y[1] = sinh( U[1] ); if (ucom1.m2 == 0) return; } d2 = Y[1]; d1 = cosh( U[1] ); if (ucom1.m1 <= 1) goto L_16; goto L_24; /* 1ST PARTIALS */ L_16: for (i = 1; i <= ucom1.n; i++) { Y[i + 1] = d1*U[i + 1]; } if (ucom1.m2 == 1) return; /* 2ND PARTIALS */ L_24: k = ucom1.n + 2; for (i = 1; i <= ucom1.n; i++) { fac = d2*U[i + 1]; for (j = 1; j <= i; j++) { Y[k] = fac*U[j + 1] + d1*U[k]; k += 1; } } return; } /* end of function */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ void /*FUNCTION*/ ducosh( double u[], double y[]) { long int i, j, k; double d1, d2, fac; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const U = &u[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* COMPUTE.. Y=COSH(U) * */ if (ucom1.m1 <= 0) { Y[1] = cosh( U[1] ); if (ucom1.m2 == 0) return; } d2 = Y[1]; d1 = sinh( U[1] ); if (ucom1.m1 <= 1) goto L_16; goto L_24; /* 1ST PARTIALS */ L_16: for (i = 1; i <= ucom1.n; i++) { Y[i + 1] = d1*U[i + 1]; } if (ucom1.m2 == 1) return; /* 2ND PARTIALS */ L_24: k = ucom1.n + 2; for (i = 1; i <= ucom1.n; i++) { fac = d2*U[i + 1]; for (j = 1; j <= i; j++) { Y[k] = fac*U[j + 1] + d1*U[k]; k += 1; } } return; } /* end of function */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ void /*FUNCTION*/ duatan( double u1[], double y1[]) { long int i, j, k; double d1, fac, fac2; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const U1 = &u1[0] - 1; double *const Y1 = &y1[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* COMPUTE.. Y= ATAN(U) * */ switch (IARITHIF(ucom1.m1 - 1)) { case -1: goto L_30; case 0: goto L_32; case 1: goto L_36; } L_30: Y1[1] = atan( U1[1] ); if (ucom1.m2 == 0) return; L_32: d1 = ONE/(ONE + SQ(U1[1])); for (i = 1; i <= ucom1.n; i++) { Y1[i + 1] = d1*U1[i + 1]; } if (ucom1.m2 == 1) return; goto L_38; L_36: d1 = ONE/(ONE + SQ(U1[1])); L_38: fac = -(U1[1] + U1[1]); k = ucom1.n + 2; for (i = 1; i <= ucom1.n; i++) { fac2 = fac*Y1[i + 1]; for (j = 1; j <= i; j++) { Y1[k] = d1*U1[k] + fac2*Y1[j + 1]; k += 1; } } return; } /* end of function */ void /*FUNCTION*/ duatn2( double v2[], double u2[], double y2[]) { long int i, j, k; double fac1, fac2, r; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const U2 = &u2[0] - 1; double *const V2 = &v2[0] - 1; double *const Y2 = &y2[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* COMPUTE Y= ATAN2(V,U) = ATAN(V/U) with 4-quadrant resolution. * */ if (ucom1.m1 == 0) { Y2[1] = atan2( V2[1], U2[1] ); if (ucom1.m2 == 0) return; } if (fabs( V2[1] ) >= fabs( U2[1] )) { r = U2[1]/V2[1]; fac2 = ONE/(V2[1] + r*U2[1]); fac1 = r*fac2; } else { r = V2[1]/U2[1]; fac1 = ONE/(V2[1]*r + U2[1]); fac2 = r*fac1; } if (ucom1.m1 == 2) goto L_54; /* FIRST PARTIALS */ for (i = 1; i <= ucom1.n; i++) { Y2[i + 1] = fac1*V2[i + 1] - fac2*U2[i + 1]; } if (ucom1.m2 == 1) return; /* SECOND PARTIALS */ L_54: k = ucom1.n + 2; for (i = 1; i <= ucom1.n; i++) { for (j = 1; j <= i; j++) { Y2[k] = fac1*(V2[k] - U2[i + 1]*Y2[j + 1] - U2[j + 1]* Y2[i + 1]) - fac2*(U2[k] + V2[i + 1]*Y2[j + 1] + V2[j + 1]* Y2[i + 1]); k += 1; } } return; } /* end of function */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ void /*FUNCTION*/ duasin( double u1[], double y1[]) { /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const U1 = &u1[0] - 1; double *const Y1 = &y1[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* COMPUTE.. Y= ASIN(U) * */ duacs( FALSE, u1, y1 ); return; } /* end of function */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ void /*FUNCTION*/ duacos( double u1[], double y1[]) { /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const U1 = &u1[0] - 1; double *const Y1 = &y1[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* COMPUTE.. Y= ACOS(U) * */ duacs( TRUE, u1, y1 ); return; } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ duacs( LOGICAL32 acosin, double u1[], double y1[]) { long int i, j, k; double fac, fac2, s1; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const U1 = &u1[0] - 1; double *const Y1 = &y1[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* COMPUTE.. Y= ACOS(U) if ACOSIN .eq. .true. and * Y= ASIN(U) if ACOSIN .eq. .false. * */ if (ucom1.m1 == 0) { if (acosin) { Y1[1] = acos( U1[1] ); } else { Y1[1] = asin( U1[1] ); } } if (ucom1.m2 == 0) return; s1 = ONE - SQ(U1[1]); if (s1 == ZERO) { /* Error condition. */ if (acosin) { ermsg( "DUACOS", 1, 0, "Deriv of ACOS(x) is infinite at x = -1 or +1" , '.' ); } else { ermsg( "DUASIN", 1, 0, "Deriv of ASIN(x) is infinite at x = -1 or +1" , '.' ); } return; } fac = ONE/sqrt( s1 ); if (acosin) fac = -fac; if (ucom1.m1 <= 1 && 1 <= ucom1.m2) { /* First partials */ for (i = 1; i <= ucom1.n; i++) { Y1[i + 1] = fac*U1[i + 1]; } } if (ucom1.m2 == 1) return; /* Second partials */ k = ucom1.n + 2; for (i = 1; i <= ucom1.n; i++) { fac2 = U1[1]*Y1[i + 1]; for (j = 1; j <= i; j++) { Y1[k] = fac*(U1[k] + fac2*Y1[j + 1]); k += 1; } } return; } /* end of function */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ void /*FUNCTION*/ dutan( double u[], double y[]) { long int i, j, k; double d1, d2, fac2; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const U = &u[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* COMPUTE.. Y=TAN(U) * */ if (ucom1.m1 <= 0) { Y[1] = tan( U[1] ); if (ucom1.m2 == 0) return; } d1 = powi(ONE/cos( U[1] ),2); d2 = TWO*Y[1]*d1; if (ucom1.m1 <= 1) goto L_216; goto L_224; /* 1ST PARTIALS */ L_216: for (i = 1; i <= ucom1.n; i++) { Y[i + 1] = d1*U[i + 1]; } if (ucom1.m2 == 1) return; /* 2ND PARTIALS */ L_224: k = ucom1.n + 2; for (i = 1; i <= ucom1.n; i++) { fac2 = d2*U[i + 1]; for (j = 1; j <= i; j++) { Y[k] = fac2*U[j + 1] + d1*U[k]; k += 1; } } return; } /* end of function */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ void /*FUNCTION*/ dutanh( double u[], double y[]) { long int i, j, k; double d1, d2, fac2; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const U = &u[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* COMPUTE.. Y=TANH(U) * */ if (ucom1.m1 <= 0) { Y[1] = tanh( U[1] ); if (ucom1.m2 == 0) return; } d1 = powi(ONE/cosh( U[1] ),2); d2 = -TWO*Y[1]*d1; if (ucom1.m1 <= 1) goto L_216; goto L_224; /* 1ST PARTIALS */ L_216: for (i = 1; i <= ucom1.n; i++) { Y[i + 1] = d1*U[i + 1]; } if (ucom1.m2 == 1) return; /* 2ND PARTIALS */ L_224: k = ucom1.n + 2; for (i = 1; i <= ucom1.n; i++) { fac2 = d2*U[i + 1]; for (j = 1; j <= i; j++) { Y[k] = fac2*U[j + 1] + d1*U[k]; k += 1; } } return; } /* end of function */