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