// RKintsubtest.cc // This file contains functions for use with RKint#test.cc #include #include #include #include #include #include #include #include #include "cpgplot.h" #include "RKtest.h" #include "RKinttest.h" // freqproc // Takes velocity amplitude data from integration routines and // writes it to a file in way that it easier for the graphing program // to process. Output forms five columns: harmonic number n, // Re[v_n], Im[v_n], |v_n|, 20*log_10(|v_n|). // // Author : Ron Kumon // Written : 11 Dec 1996 // Modified: 02 Feb 1997 for direct inclusion in RKint // 22 Feb 1997 to allow passing of the number of harmonics void freqproc(const char *outfile, const int NM, const double *V, float *normdB) { const int PRECISION=12; // Define output formatting parameters const int COLW=20; // Set column width double Vr, Vi; // Harmonics amplitudes (real and imag) double norm; // Norm of amplitude (straight) cout << "Writing Outfile= " << outfile << endl; ofstream OutputFile; // Declaration of output file object OutputFile.open(outfile, ios::out); // Open data file if (!OutputFile) { // Error checking cerr << "File could not be opened" << endl; exit(1); } OutputFile.precision(PRECISION); // Set output precision OutputFile.setf(ios::scientific); // Set scientific notation for (int n=1;n<=NM;n++) { Vr=V[2*n-1]; Vi=V[2*n]; norm=sqrt(Vr*Vr+Vi*Vi); if (norm <= 1e-40) normdB[n]=-800.0; else normdB[n]=20*log10(norm); OutputFile << setw(4) << n << setw(COLW) << Vr << setw(COLW) << Vi << setw(COLW) << norm << setw(COLW) << normdB[n] << endl; } OutputFile.close(); // Close data file } // normdBproc // Takes velocity amplitude data from integration routines // and computes normdB=20*log_10(|v_n|). // // Author : Ron Kumon // Written : 25 Sep 1997 void normdBproc(const int NM, const double *V, float *normdB) { double Vr, Vi; // Harmonics amplitudes (real and imag) double norm; // Norm of amplitude (straight) for (int n=1;n<=NM;n++) { Vr=V[2*n-1]; Vi=V[2*n]; norm=sqrt(Vr*Vr+Vi*Vi); if (norm <= 1e-40) normdB[n]=-800.0; else normdB[n]=20*log10(norm); } } // integrateFixed // This routine integrates from x1 to x2 using rkfixed.dc It also // computes the amount of time taken by the integration and // writes out general integration statistics to a file // // Author : Ron Kumon // Written : 26 Feb 1997 void integrateFixed(const char *outputfile, const char *suffix, const int nm, const int nm2, const double x1, const double x2, const double dx, double *V, double *totaltime, float const *nvec, float *normdB, float *normdBold) { const int STRLEN=81; // String length for filenames clock_t clock1, clock0; // Clock variables for timing time_t now; // Variable for current time and date double inttime; // Time for each integration char datafile[STRLEN]; // Output file + suffix for position data time(&now); // Find current time and display cout << " Starting time: " << ctime(&now); clock0=clock(); // Start time rkfixed(V, nm2, x1, x2, dx, derivs, nvec, normdB, normdBold); clock1=clock(); // End time inttime=((double)(clock1-clock0)/CLOCKS_PER_SEC); // Total time cout.precision(2); // Set precision cout.unsetf(ios::scientific); // Unset scientific notation cout.setf(ios::fixed); // Set fixed notation cout << " Time to integrate from " << x1 << " to " << x2 << " was approximately " << inttime << " s." << endl; *totaltime=*totaltime+inttime; ofstream StatFile; // Decl. of output file object StatFile.open(statfilename, ios::app); // Append to data file if (!StatFile) { // Error checking cerr << "File could not be opened" << endl; exit(1); } StatFile.precision(3); // Set precision StatFile.setf(ios::fixed); // Set fixed format StatFile << "To x= " << setw(6) << x2 << endl; StatFile.close(); // cout << endl; // for (int n=1;n<=nm2;n++) { // cout << "V[" << n << "]= " << V[n] << endl; // } strcpy(datafile,outputfile); // Copy root name to string strcat(datafile,suffix); // Add datafile suffix freqproc(datafile, nm, V, normdB); // Write data to file } void createSuffix(const char *datastr, char *newstr) // This function takes a string that contains a number, prepends a // prefix and postpends a suffix and returns the new string. It is // used to make the suffixes for the position data files. { const char prefix[]="-"; const char suffix[]=".dat"; strcpy(newstr,prefix); strcat(newstr,datastr); strcat(newstr,suffix); } void openGraph(const float xmax) { if (cpgopen("/XWIN") !=1) // Open XWindow device window exit(EXIT_FAILURE); cpgask(1); // Turn on window prompting cpgpage(); // Clear page cpgvstd(); // Define standard viewport cpgswin(1.0, xmax, -200.0, 0.0); // Define window cpglab("n", "V\\dn\\u/V\\d1\\u [dB]", "Spectrum"); // Set graph label } void eraseAndDrawGraph(const int N, const float *x, const float *y, const float t, float *yold) { static char tstringold[15]=""; // Previous parameter string value char tstring[15]; // Current parameter string value sprintf(tstring, "x = %6.3f", t); // Convert float to string // cout << "-------------------------------------------------" << endl; // cout << "Time in= " << tstring // << " Time old= " << tstringold << endl; // cout << "x[1]= " << x[1] // << " y[1]= " << y[1] // << " yold[1]= " << yold[1] << endl; cpgbbuf(); // Begin graphics buffering cpgsci(0); // Set color index to black cpgline(N, &x[1], &yold[1]); // Erase old line cpgtext(N/40+1, 10.0, tstringold); // Erase old text cpgsci(1); // Set color index to white cpgline(N, &x[1], &y[1]); // Draw new line cpgbox("bcnst", 0.0, 0, "bcnst", 0.0, 0); // Draw new frame cpgtext(N/40+1, 10.0, tstring); // Draw new text cpgebuf(); // End graphics buffering strcpy(tstringold, tstring); // Transfer current to previous for(int n=1; n<=N; n++){ yold[n] = y[n]; } // cout << "Time in= " << tstring // << " Time old= " << tstringold << endl; // cout << "x[1]= " << x[1] // << " y[1]= " << y[1] // << " yold[1]= " << yold[1] << endl; // cout << "-------------------------------------------------" << endl; } void closeGraph(void) { cpgclos(); // Close window }