/*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 "sucomp.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*/ susetn( long nin, long m1in, long m2in) { /* Base name of file: SUCOMP * Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 2001-06-18 SUCOMP Krogh Replaced "M1 .eq.. 0" with "M1 .eq. 0" *>> 1996-04-30 SUCOMP Krogh Removed redundant save statement. *>> 1994-11-02 SUCOMP Krogh Changes to use M77CON *>> 1994-08-04 SUCOMP 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 SUCOMP Lawson Initial code. * * The file SUCOMP 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 SUSETN ( N, M1, M2) * subroutine SUGETN ( N, M1, M2, L1, L2) * subroutine SUSET ( VAL, KEY, Y9) * * subroutine SUPRO ( U,V,Y) * subroutine SUQUO ( U1,V1,Y1) * subroutine SUSUM ( U3,V3,Y3) * subroutine SUDIF ( U4,V4,Y4) * subroutine SUSUM1 ( C5,V5,Y5) * subroutine SUDIF1 ( C6,V6,Y6) * subroutine SUPRO1 ( C7,V7,Y7) * subroutine SUQUO1 ( C8,V8,Y8) * * subroutine SUSQRT (U,Y) * subroutine SUEXP (U1,Y1) * subroutine SULOG (U2,Y2) * subroutine SUPWRI (MPWR,U3,Y3) * * subroutine SUSIN (U,Y) * subroutine SUCOS (U ,Y ) * subroutine SUSINH (U ,Y ) * subroutine SUCOSH (U ,Y ) * subroutine SUATAN (U1,Y1) * subroutine SUATN2 (V2,U2,Y2) * subroutine SUASIN (U, Y) * subroutine SUACOS (U, Y) * subroutine SUTAN (U, Y) * subroutine SUTANH (U, Y) * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *--S 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 SUSETN ( 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*/ sugetn( 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*/ suset( float val, long key, float y9[]) { long int i, j, k; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *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*/ supro( float u[], float v[], float y[]) { long int i, j, k; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U = &u[0] - 1; float *const V = &v[0] - 1; float *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*/ suquo( float u1[], float v1[], float y1[]) { long int i, j, k; float vinv; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U1 = &u1[0] - 1; float *const V1 = &v1[0] - 1; float *const Y1 = &y1[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* COMPUTE.. Y=U/V * */ if (V1[1] == ZERO) { ermsg( "SUQUO", 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*/ susum( float u3[], float v3[], float y3[]) { long int k; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U3 = &u3[0] - 1; float *const V3 = &v3[0] - 1; float *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*/ sudif( float u4[], float v4[], float y4[]) { long int k; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U4 = &u4[0] - 1; float *const V4 = &v4[0] - 1; float *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*/ susum1( float c5, float v5[], float y5[]) { long int k; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const V5 = &v5[0] - 1; float *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*/ sudif1( float c6, float v6[], float y6[]) { long int k; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const V6 = &v6[0] - 1; float *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*/ supro1( float c7, float v7[], float y7[]) { long int k; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const V7 = &v7[0] - 1; float *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*/ suquo1( float c8, float v8[], float y8[]) { long int i, j, k; float fac; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const V8 = &v8[0] - 1; float *const Y8 = &y8[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* COMPUTE Y=C/V WHERE C IS A CONSTANT * */ if (V8[1] == ZERO) { ermsg( "SUQUO1", 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*/ susqrt( float u[], float y[]) { long int i, j, k; float d1, fac, fac2; /* 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=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] = sqrtf( U[1] ); if (ucom1.m2 == 0) return; L_6: if (Y[1] == ZERO) { ermsg( "SUSQRT", 2, 0, "Deriv of sqrtf(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*/ suexp( float u1[], float y1[]) { long int i, j, k; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U1 = &u1[0] - 1; float *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] = expf( 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*/ sulog( float u2[], float y2[]) { long int i, j, k; float d1; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U2 = &u2[0] - 1; float *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] = logf( 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*/ supwri( long mpwr, float u3[], float y3[]) { long int i, j, k, k1, k2; float fac, fac2, fm; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U3 = &u3[0] - 1; float *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( "SUPWRI", 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] = powif(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*/ susin( float u[], float y[]) { long int i, j, k; float d1, d2, fac; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U = &u[0] - 1; float *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* computes Y=SIN(U) * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ if (ucom1.m1 <= 0) { Y[1] = sinf( U[1] ); if (ucom1.m2 == 0) return; } d2 = -Y[1]; d1 = cosf( 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*/ sucos( float u[], float y[]) { long int i, j, k; float d1, d2, fac; /* 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=COS(U) * */ if (ucom1.m1 <= 0) { Y[1] = cosf( U[1] ); if (ucom1.m2 == 0) return; } d2 = -Y[1]; d1 = -sinf( 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*/ susinh( float u[], float y[]) { long int i, j, k; float d1, d2, fac; /* 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=SINH(U) * */ if (ucom1.m1 <= 0) { Y[1] = sinhf( U[1] ); if (ucom1.m2 == 0) return; } d2 = Y[1]; d1 = coshf( 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*/ sucosh( float u[], float y[]) { long int i, j, k; float d1, d2, fac; /* 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=COSH(U) * */ if (ucom1.m1 <= 0) { Y[1] = coshf( U[1] ); if (ucom1.m2 == 0) return; } d2 = Y[1]; d1 = sinhf( 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*/ suatan( float u1[], float y1[]) { long int i, j, k; float d1, fac, fac2; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U1 = &u1[0] - 1; float *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] = atanf( 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*/ suatn2( float v2[], float u2[], float y2[]) { long int i, j, k; float fac1, fac2, r; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U2 = &u2[0] - 1; float *const V2 = &v2[0] - 1; float *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] = atan2f( V2[1], U2[1] ); if (ucom1.m2 == 0) return; } if (fabsf( V2[1] ) >= fabsf( 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*/ suasin( float u1[], float y1[]) { /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U1 = &u1[0] - 1; float *const Y1 = &y1[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* COMPUTE.. Y= ASIN(U) * */ suacs( FALSE, u1, y1 ); return; } /* end of function */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ void /*FUNCTION*/ suacos( float u1[], float y1[]) { /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U1 = &u1[0] - 1; float *const Y1 = &y1[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* COMPUTE.. Y= ACOS(U) * */ suacs( TRUE, u1, y1 ); return; } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ suacs( LOGICAL32 acosin, float u1[], float y1[]) { long int i, j, k; float fac, fac2, s1; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const U1 = &u1[0] - 1; float *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] = acosf( U1[1] ); } else { Y1[1] = asinf( U1[1] ); } } if (ucom1.m2 == 0) return; s1 = ONE - SQ(U1[1]); if (s1 == ZERO) { /* Error condition. */ if (acosin) { ermsg( "SUACOS", 1, 0, "Deriv of ACOS(x) is infinite at x = -1 or +1" , '.' ); } else { ermsg( "SUASIN", 1, 0, "Deriv of ASIN(x) is infinite at x = -1 or +1" , '.' ); } return; } fac = ONE/sqrtf( 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*/ sutan( float u[], float y[]) { long int i, j, k; float d1, d2, fac2; /* 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=TAN(U) * */ if (ucom1.m1 <= 0) { Y[1] = tanf( U[1] ); if (ucom1.m2 == 0) return; } d1 = powif(ONE/cosf( 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*/ sutanh( float u[], float y[]) { long int i, j, k; float d1, d2, fac2; /* 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=TANH(U) * */ if (ucom1.m1 <= 0) { Y[1] = tanhf( U[1] ); if (ucom1.m2 == 0) return; } d1 = powif(ONE/coshf( 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 */