/*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 <math.h>
#include "fcrt.h"
#include "ducomp.h"
#include <stdlib.h>
/* 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 */