/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:33:19 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv pf=,p_sucomp s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include <stdio.h>
#include <stdlib.h>
#include "p_sucomp.h"
/* program DRSUCOMP
*>> 1994-10-19 DRSUCOMP Krogh Changes to use M77CON
*>> 1994-08-04 DRSUCOMP CLL New subroutine: SUSETN
*>> 1992-02-17 CLL
*>> 1990-12-13 CLL Added demo of SUREV.
*>> 1987-12-09 DRSUCOMP Lawson Initial Code.
* Demo driver for the SUCOMP package, including SUREV.
* The SUCOMP package computes partial derivatives.
* SUREV does series reversion involving 1st and 2nd partial
* derivatives of N functions of N variables.
* ------------------------------------------------------------------
*--S replaces "?": DR?UCOMP, ?UCOMP, ?UREV, ?USETN, ?USET, ?UATN2
*-- & ?UPRO, ?USUM, ?USQRT, ?UQUO, ?UATAN, ?UCOS, ?USIN, ?COPY
* ------------------------------------------------------------------ */
/* PARAMETER translations */
#define LDIM (((NMAX + 2)*(NMAX + 1))/2)
#define NMAX 3
#define XVAL 0.1e0
#define YVAL 0.2e0
#define ZVAL 0.3e0
/* end of PARAMETER translations */
int main( )
{
long int _n, i, iwork[NMAX];
float cp[LDIM], ct[LDIM], phi[LDIM], r[LDIM], rcond, rsq[LDIM],
s[LDIM], sp[LDIM], ssq[LDIM], st[LDIM], t1[LDIM], t2[LDIM], t3[LDIM],
temp[LDIM], theta[LDIM], tu[NMAX][LDIM], ut[NMAX][LDIM], work[3][NMAX][NMAX],
x[LDIM], x1[LDIM], x2[LDIM], x3[LDIM], xsq[LDIM], y[LDIM], y2[LDIM],
ysq[LDIM], z[LDIM], z2[LDIM], zbys[LDIM], zsq[LDIM];
/* OFFSET Vectors w/subscript range: 1 to dimension */
float *const Cp = &cp[0] - 1;
float *const Ct = &ct[0] - 1;
long *const Iwork = &iwork[0] - 1;
float *const Phi = &phi[0] - 1;
float *const R = &r[0] - 1;
float *const Rsq = &rsq[0] - 1;
float *const S = &s[0] - 1;
float *const Sp = &sp[0] - 1;
float *const Ssq = &ssq[0] - 1;
float *const St = &st[0] - 1;
float *const T1 = &t1[0] - 1;
float *const T2 = &t2[0] - 1;
float *const T3 = &t3[0] - 1;
float *const Temp = &temp[0] - 1;
float *const Theta = &theta[0] - 1;
float *const X = &x[0] - 1;
float *const X1 = &x1[0] - 1;
float *const X2 = &x2[0] - 1;
float *const X3 = &x3[0] - 1;
float *const Xsq = &xsq[0] - 1;
float *const Y = &y[0] - 1;
float *const Y2 = &y2[0] - 1;
float *const Ysq = &ysq[0] - 1;
float *const Z = &z[0] - 1;
float *const Z2 = &z2[0] - 1;
float *const Zbys = &zbys[0] - 1;
float *const Zsq = &zsq[0] - 1;
/* end of OFFSET VECTORS */
/* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
* Set N, M1, and M2. */
susetn( NMAX, 0, 2 );
/* Initialize independent variables. */
suset( XVAL, 1, x );
suset( YVAL, 2, y );
suset( ZVAL, 3, z );
printf(" DRSUCOMP.. Demo driver for the SUCOMP package.\n "
"This demo first transforms (x,y,z) to (r,theta,phi),\n "
"and then transforms back, including computation of\n "
"1st and 2nd partial derivatives.\n");
printf(" \n VALUE D1 D2 D3 D11 D21 D22 D31 D32"
" D33\n");
printf(" \n X =");
for(_n=0L; _n < sizeof(x)/sizeof(float); _n++)
printf("%7.3f", x[_n]);
printf("\n Y =");
for(_n=0L; _n < sizeof(y)/sizeof(float); _n++)
printf("%7.3f", y[_n]);
printf("\n Z =");
for(_n=0L; _n < sizeof(z)/sizeof(float); _n++)
printf("%7.3f", z[_n]);
printf("\n");
/* Transform from (x,y,z) to (r,theta,phi) */
suatn2( y, x, phi );
supro( x, x, xsq );
supro( y, y, ysq );
supro( z, z, zsq );
susum( xsq, ysq, ssq );
susum( ssq, zsq, rsq );
susqrt( rsq, r );
susqrt( ssq, s );
suquo( z, s, zbys );
suatan( zbys, theta );
printf(" \n R =");
for(_n=0L; _n < sizeof(r)/sizeof(float); _n++)
printf("%7.3f", r[_n]);
printf("\n PHI =");
for(_n=0L; _n < sizeof(phi)/sizeof(float); _n++)
printf("%7.3f", phi[_n]);
printf("\n THETA =");
for(_n=0L; _n < sizeof(theta)/sizeof(float); _n++)
printf("%7.3f", theta[_n]);
printf("\n");
/* Transform back from (r,theta,phi) to (x,y,z) */
sucos( phi, cp );
sucos( theta, ct );
supro( cp, ct, temp );
supro( temp, r, x2 );
susin( phi, sp );
supro( sp, ct, temp );
supro( temp, r, y2 );
susin( theta, st );
supro( st, r, z2 );
printf(" \n X =");
for(_n=0L; _n < sizeof(x2)/sizeof(float); _n++)
printf("%7.3f", x2[_n]);
printf("\n Y =");
for(_n=0L; _n < sizeof(y2)/sizeof(float); _n++)
printf("%7.3f", y2[_n]);
printf("\n Z =");
for(_n=0L; _n < sizeof(z2)/sizeof(float); _n++)
printf("%7.3f", z2[_n]);
printf("\n");
/* Set data to call SUREV.
* */
scopy( LDIM, r, 1, &ut[0][0], 1 );
scopy( LDIM, phi, 1, &ut[1][0], 1 );
scopy( LDIM, theta, 1, &ut[2][0], 1 );
tu[0][0] = X[1];
tu[1][0] = Y[1];
tu[2][0] = Z[1];
surev( (float*)ut, (float*)tu, LDIM, &rcond, iwork, (float*)work );
printf(" \n To demo SUREV we store (R, PHI, THETA) including\n"
" the first and second derivatives w.r.t. (X,Y,Z) into UT(),\n"
" and set TU() = (X, Y, Z).\n Then compute (TU1,TU2,TU3) using SUREV.\n"
" For comparison compute (T1,T2,T3) using the known functional\n"
" definition of (X, Y, Z) as a function of (R, PHI, THETA).\n");
printf(" \n TU1 =");
for (i = 1; i <= LDIM; i++)
{
printf("%7.3f", tu[0][i - 1]);
}
printf("\n TU2 =");
for (i = 1; i <= LDIM; i++)
{
printf("%7.3f", tu[1][i - 1]);
}
printf("\n TU3 =");
for (i = 1; i <= LDIM; i++)
{
printf("%7.3f", tu[2][i - 1]);
}
printf("\n");
/* For comparison set X1, X2, X3, and transform them to
* T1, T2, T3. These are the same operations as
* transforming from (r,theta,phi) to (x,y,z).
* */
suset( R[1], 1, x1 );
suset( Phi[1], 2, x2 );
suset( Theta[1], 3, x3 );
sucos( x2, cp );
sucos( x3, ct );
supro( cp, ct, temp );
supro( temp, x1, t1 );
susin( x2, sp );
supro( sp, ct, temp );
supro( temp, x1, t2 );
susin( x3, st );
supro( st, x1, t3 );
printf(" \n T1 =");
for(_n=0L; _n < sizeof(t1)/sizeof(float); _n++)
printf("%7.3f", t1[_n]);
printf("\n T2 =");
for(_n=0L; _n < sizeof(t2)/sizeof(float); _n++)
printf("%7.3f", t2[_n]);
printf("\n T3 =");
for(_n=0L; _n < sizeof(t3)/sizeof(float); _n++)
printf("%7.3f", t3[_n]);
printf("\n");
exit(0);
} /* end of function */