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