/* CF5.C Program CF5.C calculates Stumpff's c-functions c(n,x) for integers n = 0, 1, 2, 3, 4, and 5, for an input argument x. It computes the c-functions by series and recursion, and then checks the calculations using the mathematically equivalent circular or hyperbolic function representations. Input is from the keyboard and output is to the screen. CF5.C was checked out via Borland's C++ Builder 5. It is based upon the program CF5.PAS, as implemented in Borland's Turbo Pascal 5.0 for IBM PCs and compatibles. Roger L. Mansfield, 6 March 2019. */ #include #include #include /* Define type "REAL" as double. */ typedef double REAL; /* Declare function and procedure prototypes. */ REAL ABS(REAL X); REAL SQR(REAL X); REAL SIN(REAL X); REAL COS(REAL X); REAL SQRT(REAL X); REAL COSH(REAL X); REAL SINH(REAL X); char CH; FILE *OUT; REAL C(int N, REAL X); void CFUNCS(REAL Y, REAL *C0, REAL *C1, REAL *C2, REAL *C3, REAL *C4, REAL *C5); int main() { /* Declare main program local variables. */ REAL X, C0, C1, C2, C3, C4, C5; OUT = fopen("cf5_out.txt","w"); if (OUT == NULL) { puts("Cannot open cf5_out.txt for output."); goto a999; } /* Prompt for and accept keyboard input. */ printf ( "C-FUNCTION VALUES FOR ARGUMENTS X\n"); printf ( "\n"); printf ( "(See file cf5_out.txt for output)\n"); fprintf(OUT,"C-FUNCTION VALUES FOR ARGUMENTS X\n"); a100: printf( " \n"); fprintf(OUT," \n"); printf( "INPUT X: "); fprintf(OUT,"INPUT X: "); scanf ( "%lf", &X); fprintf(OUT,"%12.8lf \n", X); printf(" \n"); if (X == -1000.0) goto a999; /* Calculate c-functions for input argument X. */ CFUNCS(X,&C0,&C1,&C2,&C3,&C4,&C5); /* Output C0 - C5 for input argument X. */ printf( "C0(X) = %20.10f CHECK = %20.10f\n", C0, C(0,X)); printf( "C1(X) = %20.10f CHECK = %20.10f\n", C1, C(1,X)); printf( "C2(X) = %20.10f CHECK = %20.10f\n", C2, C(2,X)); printf( "C3(X) = %20.10f CHECK = %20.10f\n", C3, C(3,X)); printf( "C4(X) = %20.10f CHECK = %20.10f\n", C4, C(4,X)); printf( "C5(X) = %20.10f CHECK = %20.10f\n", C5, C(5,X)); printf( " \n"); fprintf(OUT,"C0(X) = %20.10f CHECK = %20.10f\n", C0, C(0,X)); fprintf(OUT,"C1(X) = %20.10f CHECK = %20.10f\n", C1, C(1,X)); fprintf(OUT,"C2(X) = %20.10f CHECK = %20.10f\n", C2, C(2,X)); fprintf(OUT,"C3(X) = %20.10f CHECK = %20.10f\n", C3, C(3,X)); fprintf(OUT,"C4(X) = %20.10f CHECK = %20.10f\n", C4, C(4,X)); fprintf(OUT,"C5(X) = %20.10f CHECK = %20.10f\n", C5, C(5,X)); fprintf(OUT," \n"); printf ( "(Type c to continue, s to stop) "); CH = getch(); if (CH == 'c') { printf("%c \n", CH); goto a100; } a999: fclose(OUT); return 0; } REAL ABS(REAL X) /* ABS returns absolute value of X. */ { return fabs(X); } REAL SQR(REAL X) /* SQR raises X to 2nd power. */ { return pow(X,2.0); } REAL SIN(REAL X) /* SIN returns sine of X. */ { return sin(X); } REAL COS(REAL X) /* COS returns cosine of X. */ { return cos(X); } REAL SQRT(REAL X) /* SQRT returns square root of X. */ { return sqrt(X); } REAL COSH(REAL X) /* COSH returns the hyperbolic cosine of X. */ { return cosh(X); } REAL SINH(REAL X) /* SINH returns the hyperbolic sine of X. */ { return sinh(X); } REAL C(int N, REAL X) /* C provides the Stumpff functions c(n,x) for n = 0, 1, 2, 3, 4, or 5. Global variables: None. This procedure was modified January 8, 1992 to use series when ABS(X) <= 0.1. It was found that sin/cos and sinh/cosh representations agree with values computed by series and recursion to better than ten significant figures. Roger L. Mansfield, 31 December 1996. */ { /* Declare local variables. */ REAL c, Y; /* Elliptical path: use circular function representation. */ if (X > 0.1) { Y = SQRT(X); switch(N) { case (0) : { c = COS(Y); break; } case (1) : { c = SIN(Y)/Y; break; } case (2) : { c = (1.0 - COS(Y))/X; break; } case (3) : { c = (Y - SIN(Y))/(X*Y); break; } case (4) : { c = (COS(Y) - (1.0 - 0.5*X))/SQR(X); break; } case (5) : { c = (SIN(Y) - (Y - 0.16666666667*X*Y))/ (SQR(X)*Y); break; } default : { } } } /* Near-parabolic or parabolic path: use series representation. */ if (ABS(X) <= 0.1) { switch(N) { case (0) : { c = (1.0-X*(1.0-X*(1.0-X*(1.0-X*(1.0-X*(1.0-X/ 132.0)/ 90.0)/ 56.0)/ 30.0)/ 12.0)/2.0); break; } case (1) : { c = (1.0-X*(1.0-X*(1.0-X*(1.0-X*(1.0-X*(1.0-X/ 156.0)/110.0)/ 72.0)/ 42.0)/ 20.0)/6.0); break; } case (2) : { c = (1.0-X*(1.0-X*(1.0-X*(1.0-X*(1.0-X*(1.0-X/ 182.0)/132.0)/ 90.0)/ 56.0)/ 30.0)/ 12.0)/ 2.0; break; } case (3) : { c = (1.0-X*(1.0-X*(1.0-X*(1.0-X*(1.0-X*(1.0-X/ 210.0)/156.0)/110.0)/ 72.0)/ 42.0)/ 20.0)/ 6.0; break; } case (4) : { c = (1.0-X*(1.0-X*(1.0-X*(1.0-X*(1.0-X*(1.0-X/ 240.0)/182.0)/132.0)/ 90.0)/ 56.0)/ 30.0)/ 24.0; break; } case (5) : { c = (1.0-X*(1.0-X*(1.0-X*(1.0-X*(1.0-X*(1.0-X/ 272.0)/210.0)/156.0)/110.0)/ 72.0)/ 42.0)/120.0; break; } default : { } } } /* Hyperbolic path: use hyperbolic function representation. */ if (-X > 0.1) { Y = SQRT(-X); switch (N) { case (0) : { c = COSH(Y); break; } case (1) : { c = SINH(Y)/Y; break; } case (2) : { c = (1.0 - COSH(Y))/X; break; } case (3) : { c = (Y - SINH(Y))/(X*Y); break; } case (4) : { c = (COSH(Y) - (1.0 - 0.5*X))/SQR(X); break; } case (5) : { c = (SINH(Y) - (Y - 0.16666666667*X*Y))/ (SQR(X)*Y); break; } default : { } } } return c; } void CFUNCS(REAL Y, REAL *C0, REAL *C1, REAL *C2, REAL *C3, REAL *C4, REAL *C5) /* This procedure calculates the following c-functions Variable Stumpff's c-function C0 c(0,x) C1 c(1,x) C2 c(2,x) C3 c(3,x) C4 c(4,X) C5 c(5,X) as needed to calculate position, velocity, and state transition matrix for uniform differential correction of a cartesian state vector. Roger L. Mansfield, 31 December 1996. */ { /* Declare local variables. */ int N; REAL X, XM; XM = 0.1; /* Reduce X by quartering until it is less than parameter XM. */ X = Y; N = 0; while (ABS(X) >= XM) { X = X/4.0; N = N + 1; } /* Calculate c-functions of reduced X by series. */ *C5 = (1.0 - X*(1.0 - X*(1.0 - X*(1.0 - X*(1.0 - X* (1.0 - X/272.0)/210.0)/156.0)/110.0)/ 72.0)/ 42.0)/120.0; *C4 = (1.0 - X*(1.0 - X*(1.0 - X*(1.0 - X*(1.0 - X* (1.0 - X/240.0)/182.0)/132.0)/ 90.0)/ 56.0)/ 30.0)/ 24.0; *C3 = 0.16666666667 - X*(*C5); *C2 = 0.5 - X*(*C4); *C1 = 1.0 - X*(*C3); *C0 = 1.0 - X*(*C2); /* Calculate c-functions of original X by recursion. */ while (N != 0 ) { N = N - 1; *C5 = 0.0625*((*C2)*(*C3) + (*C4) + (*C5)); *C4 = 0.1250*((*C2)*(*C2) + (*C4) + (*C4)); *C3 = 0.2500*((*C1)*(*C2) + (*C3)); *C2 = 0.5000* (*C1)*(*C1); *C1 = (*C1)*(*C0); *C0 = 2.0000* (*C0)*(*C0) - 1.0; } return; }