// Vibration response spectrum
#include <math.h>
#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <conio.h>
#include <vector>
#include <fstream>
#include <cmath>
#include <cstring>
#include <string>
#define MAX 20000
using namespace std;
//***********************************************
void fileoper(void);
void calc(std::vector<double> level);
void vrs(void);
void interp_psd(void);
void read_two_columns(void);
void line_test(int &kflag,const std::string & str);
void library(void);
void Dirlik_core(void);
void interp(void);
void cm0(double &rms,std::vector<double> psd);
void cm1(void);
void cm2(void);
void cm4(void);
void data_input(void);
vector<double> f;
vector<double> a;
vector<double> v;
vector<double> d;
vector<double> s;
vector<double> fn;
vector<double> ai;
vector<double> aa_psd;
vector<double> rd_psd;
vector<double> N;
vector<double> SS;
vector<double> amp;
vector<double> cumu;
vector<double> peak_range;
//***********************************************
double bex[10];
double duration;
double damp;
double level[MAX];
double octave=1./24.;
double damage,df,db,ra,grms;
double rd_rms;
double EP;
double m0,m1,m2,m4;
double oa=0.;
double ov=0.;
double od=0.;
double fstart,flast;
double f_peak=0.;
double f_freq=0.;
double pi=atan2(0.,-1.);
double tpi=2.*pi;
int fnu;
int nbex;
int in[5];
int iflag=0;
int i,ii, ic, imethod;
int j;
int last;
int m, n, nn;
int ifds;
int QQ;
double ff,aa;
char amp_label[5][20];
long num;
long numBytes=300;
char string[300];
FILE *pFile[15];
char filename[40][FILENAME_MAX];
int main()
{
 cout.precision(3);
 strcpy(amp_label[1],"acceleration");
 strcpy(amp_label[2],"velocity");
 strcpy(amp_label[3],"displacement");
 in[1]=1;
 in[2]=2;
 in[3]=3;
 cout << " \n"; cout << " vrs.cpp ver 5.3, April 25, 2014 \n"; cout << " \n"; cout << " by Tom Irvine \n"; cout << " Email: tom@vibrationdata.com \n"; cout << " \n"; cout << "This program calculates the acceleration vibration response \n"; cout << "spectrum of a single-degree-of-freedom system to a \n"; cout << "power spectral density input to the base of the system. \n\n"; cout << "The power spectral density may vary arbitrarily with frequency. \n"; cout << "Natural frequency is an independent variable. Damping is fixed.\n\n"; data_input(); vrs(); cout << "\n\n Acceleration Output files \n" << endl; cout << " vrs.txt- fn(Hz) and Response(GRMS) " << endl; cout << " tsvrs.txt- fn(Hz) and Response(3-sigma) " << endl; cout << " ns_vrs.txt- fn(Hz) and Response(n-sigma) \n " << endl; cout << " where n =sqrt(2*ln(fnT)) + 0.5772/sqrt(2*ln(fnT)) \n" << endl; if(ifds==1) { cout << "\n Fatigue Damage Spectrum - fn(Hz) and Damage Index " << endl; cout << " from acceleration (peak-valley)/2 \n" << endl; for(int jk=0;jk<nbex;jk++) { fnu=jk+20; cout <<" " <<filename[fnu] << endl; } } cout << " " << endl; cout << " vrs.out - text output" << endl; cout << "\n\n Relative Displacement Output files \n" << endl; cout << " rd_vrs.txt- fn(Hz) and Response(inch RMS) " << endl; cout << " ts_rd_vrs.txt- fn(Hz) and Response(inch 3-sigma) " << endl; cout << " ns_rd_vrs.txt- fn(Hz) and Response(inch n-sigma) " << endl; cout << "\n\n n-sigma Output file: ns.grp" << endl; cout << " \n Please call the *.txtfiles into your own plotting program.\n" << endl; return 0; } void vrs(void) { double c1, c2; double inch_rms; double rho; double sum; double t, tdr; double upper=10000.; double cc; double nnf; int k; printf("\n*********Output Parameters************************ \n"); fprintf(pFile[5],"\n*********Output Parameters***************\n"); fprintf(pFile[5],"\n damping = %5.3lf \n\n",damp); fprintf(pFile[5],"\n fn Response Response"); fprintf(pFile[5],"\n (Hz) (GRMS) (3 sigma) \n"); if(fn.size()==0) { exit(1); } printf("\n calculating.."); k=0; for(long i=0; i < fn.size(); i++) { // printf(" i=%ld fn=%8.4g \n",i,fn[i]); k++; if(k==15) { printf("."); k=0; } // absolute acceleration sum=0.; aa_psd.clear(); for(j=0; j <fn.size(); j++) { rho = fn[j]/fn[i]; tdr=2.*damp*rho; c1=0.; c2=0.; c1= pow( tdr, 2. ); c2= pow( ( 1.- pow(rho,2.)), 2.); t= (1.+ c1 ) / ( c2 + c1 ); aa_psd.push_back(t*ai[j]); // printf(" fni=%8.4g fnj=%8.4g t=%8.4g tdr=%8.4g \n",fn[i],fn[j],t,tdr); // printf(" fn=%8.4g ai=%8.4g aa_psd=%8.4g\n",fn[j],ai[j],aa_psd[j]); } cm0(grms,aa_psd); // returns grms // printf("\n out of cm0: grms=%8.4g \n",grms); // getch(); if(grms>f_peak)
		{
			f_peak=grms;
			f_freq=fn[i];
		}
		if(ifds==1)
		{
 cm1();
 cm2();
 cm4();
 EP=sqrt(m4/m2);
// printf("\n go to Dirlik \n");
// printf("\n i=%ld \n",i);
 Dirlik_core();
// printf("\n back from Dirlik \n");
// printf("\n i=%ld \n",i);
 damage=0.;
 for(int kv=0;kv<nbex;kv++) { fnu=kv+20; for(long jk=0;jk<amp.size();jk++) { damage+=pow(amp[jk],bex[kv]); } fprintf(pFile[fnu],"%8.4g \t %8.4g \n",fn[i],damage); } } // printf("\n reference 4 \n"); // printf("\n i=%ld \n",i); fprintf(pFile[3],"%lf \t %lf\n",fn[i],grms); // printf("ac %lf \t %lf\n",fn[i],grms); fprintf(pFile[4],"%lf \t %lf\n",fn[i],grms*3.); // getch(); if(fn[i]>=fstart && fn[i]<=flast) { cc=sqrt(2*log(fn[i]*duration)); nnf=cc + 0.5772/cc; fprintf(pFile[8],"%lf \t %lf\n",fn[i],grms*nnf); fprintf(pFile[10],"%lf \t %lf\n",fn[i],nnf); } fprintf(pFile[5],"%10.2lf \t %10.3lf \t %10.3lf\n",fn[i], grms, grms*3.); // relative displacement rd_psd.clear(); for(j=0; j <fn.size(); j++) { c1=0.; c2=0.; c1= pow( ( pow(fn[i],2.)-pow(fn[j],2.) ),2. ); c2= pow( ( 2.*damp*fn[i]*fn[j]),2.); t= 1. / ( c2 + c1 ); if(j==0){t/=2.;} rd_psd.push_back(t*ai[j]); } cm0(rd_rms,rd_psd); inch_rms=rd_rms*386./pow(tpi,2.); fprintf(pFile[7],"%lf \t %12.3e\n",fn[i],inch_rms); // printf("rd %lf \t %12.3e\n",fn[i],inch_rms); fprintf(pFile[11],"%lf \t %12.3e\n",fn[i],3.*inch_rms); if(fn[i]>=fstart && fn[i]<=flast) { cc=sqrt(2*log(fn[i]*duration)); nnf=cc + 0.5772/cc; fprintf(pFile[9],"%lf \t %12.3e\n",fn[i],inch_rms*nnf); } } fclose( pFile[3] ); fclose( pFile[4] ); fclose( pFile[5] ); fclose( pFile[7] ); fclose( pFile[8] ); fclose( pFile[9] ); fclose( pFile[10] ); fclose( pFile[11] ); if(ifds==1) { for(int jk=0;jk<nbex;jk++) { fnu=jk+20; fclose( pFile[fnu] ); } } printf("\n\n peak VRS at %8.4g Hz %8.4g GRMS \n",f_freq,f_peak); } void calc(std::vector<double> level)
{
		 ra=0.;
 grms=0.;
		 s.clear();
		 for(i=0; i < m; i++) { if(level[i]<1.0e-90) { printf("* f[%ld]=%8.4g level[%ld]=%10.4e \n",i,f[i],i,level[i]); printf("\n error 1: %s \n",amp_label[ii]); exit(1); } } // calculate slopes nn=m-1; double lower_f=0.0001; for(i=0; i < nn; i++) { if(f[i] < lower_f){f[i]=lower_f;} if(f[i]>f[i+1])
			 {
					printf("\n error 2: %s\n",amp_label[ii]);
					exit(1);
			 }
			 if(f[i]< lower_f) { printf("\n error 3: %s \n",amp_label[ii]); exit(1); } if(f[i+1]< lower_f) { printf("\n error 4: %s \n",amp_label[ii]); exit(1); } if(level[i]<1.0e-90) { printf("***** f[%ld]=%8.4g level[%ld]=%10.4e \n",i,f[i],i,level[i]); printf("\n error 5: %s \n",amp_label[ii]); exit(1); } s.push_back( log10( level[i+1]/ level[i] )/log10( f[i+1]/f[i] ) ); } for( i=0; i < nn; i++) { if( fabs( s[i]+1.)>1.0e-03 )
			 {
 ra+= ( level[i+1] * f[i+1]- level[i]*f[i])/( s[i]+1.);
			 }
		 else
			 {
 ra+= level[i]*f[i]*log( f[i+1]/f[i]);
			 }
		 }
 grms=sqrt(ra);
}
void fileoper(void)
{
 strcpy(filename[1], "a.psd");
 pFile[1]=fopen(filename[1], "w");
	 strcpy(filename[2], "a.int");
 pFile[2]=fopen(filename[2], "w");
	 strcpy(filename[3], "vrs.grp");
 pFile[3]=fopen(filename[3], "w");
	 strcpy(filename[4], "tsvrs.grp");
 pFile[4]=fopen(filename[4], "w");
	 strcpy(filename[5], "vrs.out");
 pFile[5]=fopen(filename[5], "w");
	 strcpy(filename[6], "vrs.in");
 pFile[6]=fopen(filename[6], "w");
	 strcpy(filename[7], "rd_vrs.grp");
 pFile[7]=fopen(filename[7], "w");
	 strcpy(filename[8], "ns_vrs.grp");
 pFile[8]=fopen(filename[8], "w");
	 strcpy(filename[9], "ns_rd_vrs.grp");
 pFile[9]=fopen(filename[9], "w");
	 strcpy(filename[10], "ns.grp");
 pFile[10]=fopen(filename[10], "w");
	 strcpy(filename[11], "ts_rd_vrs.grp");
 pFile[11]=fopen(filename[11], "w");
}
void library(void)
{
 int ilibrary=0;
 int Lflag=0;
 while(Lflag==0)
 {
 printf("\n Select PSD: \n");
 printf("\n 1=MIL-STD-1540B QTP ");
 printf("\n 2=MIL-STD-1540B ATP ");
 printf("\n 3=NAVMAT P9495 ");
 printf("\n 4=NASA CEQATR AVT \n");
 scanf("%d", &ilibrary);
 if(ilibrary==1 &#124;&#124; ilibrary==2 )
 {
 f.push_back(20.);
 a.push_back(0.0053);
 f.push_back(150.);
 a.push_back(0.04);
 f.push_back(600.);
 a.push_back(0.04);
 f.push_back(2000.);
 a.push_back(0.0036);
 Lflag=1;
 }
 if(ilibrary==1)
 {
 a[0]*=4.;
 a[1]*=4.;
 a[2]*=4.;
 a[3]*=4.;
 }
 if(ilibrary==3)
 {
 f.push_back(20.);
 a.push_back(0.01);
 f.push_back(80.);
 a.push_back(0.04);
 f.push_back(350.);
 a.push_back(0.04);
 f.push_back(2000.);
 a.push_back(0.007);
 Lflag=1;
 }
 if(ilibrary==4)
 {
 f.push_back(20.);
 a.push_back(0.01);
 f.push_back(80.);
 a.push_back(0.04);
 f.push_back(800.);
 a.push_back(0.04);
 f.push_back(2000.);
 a.push_back(0.01);
 Lflag=1;
 }
 n=f.size();
 }
 i=n;
}
void read_two_columns(void)
{
 char string_a[300];
 long i=0;
 double ff,aa;
 while( fgets(string_a,numBytes,pFile[0])> 0 )
	{
//	 printf("%s",string_a);
 int kflag=0;
 line_test(kflag,string_a);
 if(kflag==0)
 {
 if(strrchr(string_a,',' ) != NULL )
 {
 sscanf(string_a,"%lf, %lf",&ff, &aa);
 }
 else
 {
 sscanf(string_a,"%lf %lf", &ff, &aa);
 }
 f.push_back(ff);
 a.push_back(aa);
// printf(" %8.4g %8.4g \n",f[i],a[i]);
 i++;
 }
	}
	n=i;
}
void line_test(int &kflag,const std::string & str)
{
// detect header lines
		if (std::string::npos != str.find_first_of("Ee"))
		{
		 std::size_t f1=str.find("e+");
		 std::size_t f2=str.find("e-");
		 std::size_t f3=str.find("E+");
		 std::size_t f4=str.find("E-");
		 if(f1==std::string::npos && f2==std::string::npos && f3==std::string::npos && f4==std::string::npos)
		 {
		 kflag=1;
		 }
		}
 if (std::string::npos != str.find_first_of("ABCDFGHIJKLMNOPQRSTUVWXYZabcdfghijklmnopqrstuvwxyz"))
 {
 kflag=1;
 }
 if (std::string::npos != str.find_first_of("!@#$%^&*/<strong>?()=;"))
 {
 kflag=1;
 }
 if (std::string::npos == str.find_first_of("0123456789"))
 {
 kflag=1;
 }
}
void Dirlik_core(void)
{
// printf("\n in core \n");
 double gamma;
 double D1,D12,D2,D3;
 double Q,R,R2;
 double maxS;
 double ds;
 double x;
 double Z;
 x=(m1/m0)*sqrt(m2/m4);
 gamma=m2/(sqrt(m0*m4));
 D1=2*(x-pow(gamma,2))/(1+pow(gamma,2));
 D12=pow(D1,2);
 R=(gamma-x-D12)/(1-gamma-D1+D12);
 R2=pow(R,2);
 D2=(1-gamma-D1+D12)/(1-R);
 D3=1-D1-D2;
 Q=1.25*(gamma-D3-D2*R)/D1;
 maxS=8*grms;
// printf(" rms %8.4g ",grms);
 ds=maxS/400.;
 int n=int(0.5+(maxS/ds));
// printf(" ds %8.4g ",ds);
// printf(" maxS %8.4g ",maxS);
// printf("\n n %ld \n",n);
 N.clear();
 SS.clear();
 double a,b;
 double t1,t2,t3;
 double p,pn,pd;
// printf(" grms=%8.4g \n",grms);
// printf(" m0=%8.4g m1=%8.4g m2=%8.4g m3=%8.4g \n",m0,m1,m2,m4);
// printf(" EP=%8.4g \n",EP);
// printf(" Q=%8.4g \n",Q);
// printf(" ds=%8.4g \n",ds);
 for(long i=0;i<n;i++) // keep long { // printf(" i=%ld \n",i); SS.push_back(i*ds); Z=SS[i]/(2*sqrt(m0)); t1=(D1/Q)*exp(-Z/Q); a=-pow(Z,2); b=2*pow(R,2); t2=(D2*Z/R2)*exp(a/b); t3=D3*Z*exp(-pow(Z,2)/2); pn=t1+t2+t3; pd=2*sqrt(m0); p=pn/pd; N.push_back(p*EP*duration); // printf(" N=%8.4g p=%8.4g EP=%8.4g duration=%8.4g \n",N[i],p,EP,duration); } cumu.clear(); double area=0.; // printf("\n N.size %d \n",N.size()); for(long i=0;i<n.size();i++) { area+=N[i]*ds; cumu.push_back(area); // printf(" %ld cumu=%10.45 N=%8.4g ds=%8.4g\n",i,cumu[i],N[i],ds); } // print(' i=%d N=%g area=%g ds=%g ' %(i,N[i],area,ds)) // printf("\n cumu=%8.4g \n",cumu.back()); num=long(0.5+cumu.back()); // printf("\n num %ld \n",num); // vq1 = interp(xq,cumu,S) // print(' len vq1 %d ' %(len(vq1))) // peak_range=array(vq1) // peak_range=sorted(peak_range, reverse=True) // nn=len(peak_range) interp(); amp.clear(); // printf("\n pr size %d \n",peak_range.size()); double maxamp=0.; for(long i=0;i<peak_range.size();i++) { amp.push_back(peak_range[i]/2.); if(amp[i]>maxamp)
 {
 maxamp=amp[i];
 }
 }
// printf("\n amp size %d max amp=%8.4g maxS=%8.4g\n",amp.size(),maxamp,maxS);
// printf(" end DDD");
}
void interp(void)
{
 double x,l;
 double c1,c2;
 double ttt;
 peak_range.clear();
	peak_range.push_back(0.);
 long jfirst=1;
 for(long i=0;i<num;i++) // keep long { peak_range.push_back(0.); } for(long i=1; i<num; i++) { ttt = i; // printf(" ttt=%8.4g cumu size %ld\n",ttt,cumu.size()); for(long j=jfirst; j<cumu.size(); j++) { if( ttt>= cumu[j-1] && ttt <= cumu[j]) { x=ttt-cumu[j-1]; l=cumu[j]-cumu[j-1]; c2=x/l; c1=1.-c2; peak_range[i]=(c1*SS[j-1] + c2*SS[j]); // printf(" pr=%8.4g \n",peak_range[i]); jfirst=j-1; if(jfirst<0) { jfirst=0; } break; } } } } void cm0(double &rms,std::vector<double> psd)
{
// m0=m0+aa_psd*df;
 double ra=0.;
 double s;
// printf(" (psd.size()-1)=%ld \n",(psd.size()-1));
 for(long i=0; i < (psd.size()-1); i++) { // printf(" i=%ld fn=%8.4g fstart=%8.4g flast=%8.4g \n",i,fn[i],fstart,flast); if(fn[i]>=fstart && fn[i+1]<=flast) { s=( log( psd[i+1]/psd[i] )/log( fn[i+1]/fn[i] ) ); if(s < -1.0001 &#124;&#124; s> -0.9999 )
 {
 ra+= ( psd[i+1] * fn[i+1]- psd[i]*fn[i])/( s+1.);
 }
 else
 {
 ra+= psd[i]*fn[i]*log( fn[i+1]/fn[i]);
 }
 }
 }
 m0=ra;
 rms=sqrt(ra);
}
void cm1(void)
{
// m1=m1+aa_psd*fn[j]*df;
 double ra=0.;
 double s;
 vector<double> amp;
 for(long i=0; i < aa_psd.size(); i++) { amp.push_back( aa_psd[i]*fn[i] ); } for(long i=0; i < (aa_psd.size()-1); i++) { if(fn[i]>=fstart && fn[i+1]<=flast) { s=( log( amp[i+1]/amp[i] )/log( fn[i+1]/fn[i] ) ); if(s < -1.0001 &#124;&#124; s> -0.9999 )
 {
 ra+= ( amp[i+1] * fn[i+1]- amp[i]*fn[i])/( s+1.);
 }
 else
 {
 ra+= amp[i]*fn[i]*log( fn[i+1]/fn[i]);
 }
 }
 }
 m1=ra;
}
void cm2(void)
{
// m2=m2+aa_psd*pow(fn[j],2)*df;
 double ra=0.;
 double s;
 vector<double> amp;
 for(long i=0; i < aa_psd.size(); i++) { amp.push_back( aa_psd[i]*pow(fn[i],2.) ); } for(long i=0; i < (aa_psd.size()-1); i++) { if(fn[i]>=fstart && fn[i+1]<=flast) { s=( log( amp[i+1]/amp[i] )/log( fn[i+1]/fn[i] ) ); if(s < -1.0001 &#124;&#124; s> -0.9999 )
 {
 ra+= ( amp[i+1] * fn[i+1]- amp[i]*fn[i])/( s+1.);
 }
 else
 {
 ra+= amp[i]*fn[i]*log( fn[i+1]/fn[i]);
 }
 }
 }
 m2=ra;
}
void cm4(void)
{
// m4=m4+aa_psd*pow(fn[j],4)*df;
 double ra=0.;
 double s;
 vector<double> amp;
 for(long i=0; i < aa_psd.size(); i++) { amp.push_back( aa_psd[i]*pow(fn[i],4.) ); } for(long i=0; i < (aa_psd.size()-1); i++) { if(fn[i]>=fstart && fn[i+1]<=flast) { s=( log( amp[i+1]/amp[i] )/log( fn[i+1]/fn[i] ) ); if(s < -1.0001 &#124;&#124; s> -0.9999 )
 {
 ra+= ( amp[i+1] * fn[i+1]- amp[i]*fn[i])/( s+1.);
 }
 else
 {
 ra+= amp[i]*fn[i]*log( fn[i+1]/fn[i]);
 }
 }
 }
 m4=ra;
}
//***************************************************************
void interp_psd(void)
{
 long nbreak=a.size();
 vector<double> slope;
	ai.push_back(a[0]);
	for(i=0; i< (nbreak-1);i++) { slope.push_back(log(a[i+1]/a[i])/log(f[i+1]/f[i])); // printf(" i=%ld slope=%8.4g \n",i,slope[i]); } // printf("\n fn.size()=%ld \n\n",fn.size()); // printf(" fn=%8.4g ai=%8.4g \n",fn[0],ai[0]); for(i=1; i < fn.size(); i++) { int jflag=0; for(j=0; j < (nbreak-1); j++) { if( ( fn[i]>= f[j] ) && ( fn[i] <= f[j+1] ) ) { ai.push_back(a[j]*pow( ( fn[i] / f[j] ), slope[j] )); // printf(" fn=%8.4g ai=%8.4g \n",fn[i],ai[i]); jflag=1; break; } } if(jflag==0) { ai.push_back(0.); // printf(" fn=%8.4g ai=%8.4g \n",fn[i],ai[i]); } } } void data_input(void) { cout << " Select data input method \n"; cout << " 1=keyboard 2=file 3=library PSD \n"; cin>> imethod;
 if(imethod==1)
 {
 cout << " How Many Breakpoints (minimum=2)? \n"; cin>> n;
 if(n < 2) { cout << "Input error \n"; } cout << " Enter Points (free-format) "<< endl; cout << " Freq(Hz) Level(G^2/Hz) "<< endl; for( i=0; i<n; i++) { scanf("%lf %lf",&ff,&aa); f.push_back(ff); a.push_back(aa); } } if(imethod==2) { printf( " The file must contain two columns: \n"); printf( " Freq(Hz) Level(G^2/Hz) \n"); printf( " \n An error message will occur if the file does not exist.\n"); printf( " Input filename \n"); scanf("%s",filename[0]); pFile[0] = fopen(filename[0], "rb"); while(pFile[0] == NULL ) { printf("\n Failed to open file: %s \n", filename[0]); printf("\n Please enter the input filename: \n"); scanf("%s",filename[0]); pFile[0] = fopen(filename[0], "rb"); } printf(" File: %s opened. \n", filename[0]); read_two_columns(); fclose( pFile[0]); } if(imethod==3) { library(); } printf("\n"); printf("\n Enter duration (sec) "); scanf("%lf",&duration); fileoper(); m=f.size(); for(i=0; i<m; i++) { if((i>0) && (f[i] < f[i-1])) { printf(" Error: frequencies must ascend. \n"); exit(1); } if(f[i] < 0.) { printf(" all input numbers must be> 0. \n");
 printf(" %g %g \n",f[i],a[i]);
 exit(1);
 }
 if(a[i] < 0.) { printf(" all input numbers must be> 0. \n");
 printf(" %g %g \n",f[i],a[i]);
 exit(1);
 }
 fprintf( pFile[1],"%lf \t %lf \n",f[i],a[i]);
 fprintf( pFile[6],"%12.4g \t %12.4g \n",f[i],a[i]);
 }
 fclose(pFile[1]);
 fclose(pFile[6]);
 fstart=f[0];
 flast=f.back();
	double f1;
//	printf("\n\n fstart=%8.4g \n\n",fstart);
	f1=fstart;
	if(f1>5)
	{
	 f1=5;
	}
	fn.push_back(f1);
	i=1;
	while(1)
	{
	 ff=fn[i-1]*pow(2.,octave);
	 if(ff>=flast)
	 {
	 fn.push_back(flast);
	 break;
	 }
	 else
	 {
	 fn.push_back(ff);
	 }
	 i++;
	}
// printf("\n fn.size()= %ld \n",fn.size());
//*******************************************************************************
	long idamp;
	printf("\n Enter damping format. \n 1= damping ratio 2= Q \n");
	scanf("%ld",&idamp);
	if(idamp == 1 )
	{
	 printf("\n Enter damping ratio (typically 0.05) ");
 scanf("%lf",&damp);
 QQ=int(1./(2.*damp));
	}
	else
 {
	 printf("\n Enter amplification factor Q (typically 10) ");
 scanf("%lf",&damp);
 QQ=int(damp);
	 damp = 1./(2.*damp);
	}
//*******************************************************************************
 printf("\n Calculate fatigue damage spectrum? 1=yes 2=no \n");
 scanf("%d",&ifds);
 int fnu;
 if(ifds==1)
 {
 int ns;
 printf("\n Enter the number of fatigue exponents \n");
 scanf("%d",&nbex);
 if(nbex>10)
 {
 nbex=10;
 }
 for(int jk=0; jk<nbex; jk++) { fnu=jk+20; printf("\n Enter fatigue exponent %d\n",jk+1); scanf("%lf",&bex[jk]); if(imethod==2) { i = 0; ns = 0; while(filename[0][i] != '0円') // get length of filename { i++; } for(ii = 0; i> -1; i--)
 {
 // look for extension working backwards
 if(filename[0][i] == '.')
 {
 ns = i; // char # of exension
 break;
 }
 }
 memcpy(filename[fnu], filename[0], ns);
 strcat(filename[fnu],"_fds_Q");
 }
 else
 {
 strcpy(filename[fnu], "fds_Q" );
 }
 char temp[40];
 sprintf(temp,"%d_b%g.txt",QQ,bex[jk]);
 strcat(filename[fnu], temp );
 pFile[fnu]=fopen(filename[fnu], "w");
 }
 }
// printf("\n reference 1 \n");
 ii=1;
 calc(a);
// printf("\n reference 2 \n");
 interp_psd();
 oa=grms;
 printf("\n*********Input Parameters************************ \n");
 printf(" Freq PSD Slope Slope \n");
 printf(" (Hz) (G^2/Hz) log-log (dB/octave) \n");
 printf(" %6.1lf %lf \n", f[0],a[0]);
 fprintf(pFile[5],"\n*********Input Parameters************************ \n");
 fprintf(pFile[5]," Freq PSD Slope Slope \n");
 fprintf(pFile[5]," (Hz) (G^2/Hz) log-log (dB/octave) \n");
 fprintf(pFile[5]," %6.1lf %lf \n", f[0],a[0]);
 for( i=0; i < m-1; i++) { db=( 10.*log10(2.) )*s[i]; printf(" %6.1lf %lf %lf %lf \n",f[i+1],a[i+1],s[i],db); fprintf(pFile[5]," %6.1lf %lf %lf %lf \n",f[i+1],a[i+1],s[i],db); flast=f[i+1]; } fstart=f[0]; for(i=0; i<m; i++) { if( f[i]>= 1.) // highpass filter
 {
 v.push_back((pow( 386.,2.) )*a[i]/(pow( ( tpi*f[i] ),2.)));
 }
 else
 {
 v.push_back(1.0e-60);
 }
 }
 ii=2;
 calc(v);
 ov=grms;
 for(i=0; i<m; i++) { if( f[i]>= 1.) // highpass filter
 {
 d.push_back((pow( 386.,2.) )*a[i]/(pow( ( tpi*f[i] ),4.)));
 }
 else
 {
 d.push_back(1.0e-60);
 }
 if(d[i]<1.0e-80)
 {
 d[i]=1.0e-80;
 printf("\n warning: f[%ld]=%8.4g a[%ld]=%8.4g d[%ld]=%8.4g \n",i,f[i],i,a[i],i,d[i]);
 printf("\n press any key to continue \n");
 cin.get();
 }
 }
 for(i=0; i<m; i++)
 {
 if(d[i]<1.0e-80)
 {
 printf("pre-calc f[%ld]=%8.4g a[%ld]=%8.4g d[%ld]=%8.4g \n",i,f[i],i,a[i],i,d[i]);
 exit(1);
 }
 }
 ii=3;
 calc(d);
 od=grms;
 fprintf(pFile[5],"\n Overall acceleration = %10.3lf GRMS\n", oa);
 cout << "\n Overall Acceleration = " << oa << " GRMS " << endl;
 cout << " Overall Velocity = " << ov << " (inch RMS)/sec" << endl;
 cout << " Overall Displacement = " << od << " (inch RMS)" << endl;
 cout << " Overall Displacement = " << od*3 << " (inch 3-sigma)" << endl;
 cout << "\n Input points written to file: vrs.in "<< endl;
}
</div><div class="naked_ctrl">
<form action="/index.cgi/larger-text" method="get" name="gate">
<p><a href="http://altstyle.alfasado.net">AltStyle</a> によって変換されたページ <a href="http://www.vibrationdata.com/tutorials2/vrs.cpp">(-&gt;オリジナル)</a>
/ <label>アドレス: <input type="text" name="naked_post_url" value="http://www.vibrationdata.com/tutorials2/vrs.cpp" size="22" /></label> <label>モード: <select name="naked_post_mode">
<option value="default">デフォルト</option>
<option value="speech">音声ブラウザ</option>
<option value="ruby">ルビ付き</option>
<option value="contrast">配色反転</option>
<option value="larger-text" selected="selected">文字拡大</option>
<option value="mobile">モバイル</option>
</select>
<input type="submit" value="表示" />
</p>
</form>
</div>