/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:33:13 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv pf=,p_ducomp s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include #include #include "p_ducomp.h" /* program DRDUCOMP *>> 1994-10-19 DRDUCOMP Krogh Changes to use M77CON *>> 1994-08-04 DRDUCOMP CLL New subroutine: DUSETN *>> 1992-02-17 CLL *>> 1990-12-13 CLL Added demo of DUREV. *>> 1987-12-09 DRDUCOMP Lawson Initial Code. * Demo driver for the DUCOMP package, including DUREV. * The DUCOMP package computes partial derivatives. * DUREV does series reversion involving 1st and 2nd partial * derivatives of N functions of N variables. * ------------------------------------------------------------------ *--D 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]; double 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 */ double *const Cp = &cp[0] - 1; double *const Ct = &ct[0] - 1; long *const Iwork = &iwork[0] - 1; double *const Phi = &phi[0] - 1; double *const R = &r[0] - 1; double *const Rsq = &rsq[0] - 1; double *const S = &s[0] - 1; double *const Sp = &sp[0] - 1; double *const Ssq = &ssq[0] - 1; double *const St = &st[0] - 1; double *const T1 = &t1[0] - 1; double *const T2 = &t2[0] - 1; double *const T3 = &t3[0] - 1; double *const Temp = &temp[0] - 1; double *const Theta = &theta[0] - 1; double *const X = &x[0] - 1; double *const X1 = &x1[0] - 1; double *const X2 = &x2[0] - 1; double *const X3 = &x3[0] - 1; double *const Xsq = &xsq[0] - 1; double *const Y = &y[0] - 1; double *const Y2 = &y2[0] - 1; double *const Ysq = &ysq[0] - 1; double *const Z = &z[0] - 1; double *const Z2 = &z2[0] - 1; double *const Zbys = &zbys[0] - 1; double *const Zsq = &zsq[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Set N, M1, and M2. */ dusetn( NMAX, 0, 2 ); /* Initialize independent variables. */ duset( XVAL, 1, x ); duset( YVAL, 2, y ); duset( ZVAL, 3, z ); printf(" DRDUCOMP.. Demo driver for the DUCOMP 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(double); _n++) printf("%7.3f", x[_n]); printf("\n Y ="); for(_n=0L; _n < sizeof(y)/sizeof(double); _n++) printf("%7.3f", y[_n]); printf("\n Z ="); for(_n=0L; _n < sizeof(z)/sizeof(double); _n++) printf("%7.3f", z[_n]); printf("\n"); /* Transform from (x,y,z) to (r,theta,phi) */ duatn2( y, x, phi ); dupro( x, x, xsq ); dupro( y, y, ysq ); dupro( z, z, zsq ); dusum( xsq, ysq, ssq ); dusum( ssq, zsq, rsq ); dusqrt( rsq, r ); dusqrt( ssq, s ); duquo( z, s, zbys ); duatan( zbys, theta ); printf(" \n R ="); for(_n=0L; _n < sizeof(r)/sizeof(double); _n++) printf("%7.3f", r[_n]); printf("\n PHI ="); for(_n=0L; _n < sizeof(phi)/sizeof(double); _n++) printf("%7.3f", phi[_n]); printf("\n THETA ="); for(_n=0L; _n < sizeof(theta)/sizeof(double); _n++) printf("%7.3f", theta[_n]); printf("\n"); /* Transform back from (r,theta,phi) to (x,y,z) */ ducos( phi, cp ); ducos( theta, ct ); dupro( cp, ct, temp ); dupro( temp, r, x2 ); dusin( phi, sp ); dupro( sp, ct, temp ); dupro( temp, r, y2 ); dusin( theta, st ); dupro( st, r, z2 ); printf(" \n X ="); for(_n=0L; _n < sizeof(x2)/sizeof(double); _n++) printf("%7.3f", x2[_n]); printf("\n Y ="); for(_n=0L; _n < sizeof(y2)/sizeof(double); _n++) printf("%7.3f", y2[_n]); printf("\n Z ="); for(_n=0L; _n < sizeof(z2)/sizeof(double); _n++) printf("%7.3f", z2[_n]); printf("\n"); /* Set data to call DUREV. * */ dcopy( LDIM, r, 1, &ut[0][0], 1 ); dcopy( LDIM, phi, 1, &ut[1][0], 1 ); dcopy( LDIM, theta, 1, &ut[2][0], 1 ); tu[0][0] = X[1]; tu[1][0] = Y[1]; tu[2][0] = Z[1]; durev( (double*)ut, (double*)tu, LDIM, &rcond, iwork, (double*)work ); printf(" \n To demo DUREV 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 DUREV.\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). * */ duset( R[1], 1, x1 ); duset( Phi[1], 2, x2 ); duset( Theta[1], 3, x3 ); ducos( x2, cp ); ducos( x3, ct ); dupro( cp, ct, temp ); dupro( temp, x1, t1 ); dusin( x2, sp ); dupro( sp, ct, temp ); dupro( temp, x1, t2 ); dusin( x3, st ); dupro( st, x1, t3 ); printf(" \n T1 ="); for(_n=0L; _n < sizeof(t1)/sizeof(double); _n++) printf("%7.3f", t1[_n]); printf("\n T2 ="); for(_n=0L; _n < sizeof(t2)/sizeof(double); _n++) printf("%7.3f", t2[_n]); printf("\n T3 ="); for(_n=0L; _n < sizeof(t3)/sizeof(double); _n++) printf("%7.3f", t3[_n]); printf("\n"); exit(0); } /* end of function */