#include <cstdlib>
#include <iostream>
#include <vector>
#include <math.h>
#include <stdio.h>
#include <string>
#include <string.h>
#include <stdlib.h>
#include <conio.h>
#define MAX 40000000
#define ABS_SLOPE \
double sss=0.;\
double fr;\
for(i=0; i< (nbreak-1); i++)\ {\ fr=f_sam[i+1]/f_sam[i];\ sss=0.;\ sss=log(apsd_sam[i+1]/apsd_sam[i])/log(fr);\ if(sss> slopec)\
{\
apsd_sam[i+1]=apsd_sam[i]*pow(fr,slopec);\
}\
if(sss < -slopec)\ {\ apsd_sam[i+1]=apsd_sam[i]*pow(fr,-slopec);\ }\ }\ if(initial==1 && apsd_sam[0]> apsd_sam[1])\
{\
apsd_sam[0]=apsd_sam[1];\
}\
if(final==1 && apsd_sam[nbreak-1]> apsd_sam[nbreak-2])\
{\
apsd_sam[nbreak-1]=apsd_sam[nbreak-2];\
}\
#define FSORT \
for(i=1;i<=(nbreak-2);i++)\ {\ for(j=i+1; j<=(nbreak-1); j++)\ {\ if(f_sam[j] < f_sam[i] )\ {\ temp=f_sam[j];\ f_sam[j]=f_sam[i];\ f_sam[i]=temp;\ }\ }\ }\ #define FSPACE \ int ialarm = 0;\ double fnum = (fn[n_ref-1]-fn[0])/200.;\ for(i=0; i<= (nbreak-2); i++)\ {\ if ( fabs(f_sam[i+1] - f_sam[i] ) < fnum )\ {\ ialarm = 1;\ break;\ }\ }\ if(ialarm == 1 )\ {\ double df=(fn[n_ref-1]-fn[0])/(n_ref-1);\ for(i=1; i<=(nbreak-2); i++)\ {\ f_sam[i]=f_sam[i-1]+df;\ }\ }\ f_sam[nbreak-1]=fn[n_ref-1];\ FSORT\ // Dynamic Memory Allocation for 2D Arrays #define ZERO(Q,nrows,ncols) \ for(long i=0; i<nrows; i++)\ { for(long j=0; j<ncols; j++)\ {Q[i][j]=0.;}} #define DYNAMIC_MATRIX(Q,nrows,ncols) \ double **Q;\ Q =new double* [nrows];\ for(long row=0;row<nrows;row++) {Q[row]=new double[ncols]; }\ ZERO(Q,nrows,ncols) #define DELETE_MATRIX(aaa,nrows)\ for (long i = 0; i < nrows; i++)\ {delete(aaa[i]);}\ delete(aaa); \ #define DYNAMIC_VECTOR(Q,nrows)\ double* Q = NULL;\ Q = new double[nrows]; #define EXIT \ printf("\n\n Press any key to exit.\n"); \ getch(); \ printf("\n"); \ exit(1); \ #include <time.h>
#include <sys/types.h>
#include <sys/timeb.h>
using namespace std;
class FDS
{
 private:
 double alpha;
 double alpha2;
 double E,K,C,S,Sp;
 double a1,a2;
 double b1,b2,b3;
 double cosd;
 double domegadt;
 double damage;
 double flast;
 long idamp;
 long last;
 long initial;
 long i,j,jnum;
 long index;
 long ic;
 long mh;
 long nin;
 long nnn;
 long num;
 long np;
 long n_ref;
 long numBytes;
 int algor;
 int ira,iresp;
 int iunit;
 double f1,f2;
 double ocn;
 double octave;
 double RMS;
 double sind;
 double t, t1, tb;
 double x;
 double y;
 long goal;
 long final;
 double ra;
 double s;
 double slopec;
 double record;
 double pi;
 double tpi;
 char string_a[300];
 public:
 vector<vector<float>> dr; // need space>>
 double oct;
 double octh;
 double dB;
 double m0,m1,m2,m4;
 double df1,df2,df3;
 double vf1,vf2,vf3;
 double af1,af2,af3;
 double ref;
 double scale;
 double pressure_rms;
 double pressure_rmslow;
 double vrms;
 double vrmslow;
 double drms;
 double drmslow;
 double ft_rms;
 double drmsp,vrmsp,pressure_rmsp;
 double delta;
 double duration;
 double apsd_sam_last;
 double EP;
 double dt;
 double sr;
 double xf[20];
 double xapsd[20];
 double apsd_sam[20];
 double xslope[20];
 double xffine[2000];
 double xfds[2000];
 double fatigue_damage[6][6][600];
 double reference_fds[6][6][600];
 double ratio_trial_max;
 double ratio_record;
 double fc[34];
 double fband[34];
 double spl[34];
 double sum_pressure[34];
 double octbw;
 int fcases;
 int io;
 int iqcases;
 int iy;
 int iq;
 int iopt;
 int iu;
 int fc_num;
 long nn;
 long ntrials;
 long number_lines;
 long nbreak;
 long ik;
 long ikbest;
 char temp_Q[10];
 FILE *pFile[6];
 char filename[6][FILENAME_MAX];
 char input_fds_filename[10][FILENAME_MAX];
 char output_fds_filename[10][FILENAME_MAX];
 vector<int> xq;
 vector<float> yin;
 vector<double> cumu;
 vector<double> amp;
 vector<double> peak_range;
 vector<double> f_sam;
 vector<double> fin;
 vector<double> fn;
 vector<double> xabs;
 vector<double> psdin;
 vector<double> inslope;
 vector<double> interp_psdin;
 vector<double> N;
 vector<double> SS;
	 vector<double> ppsd_samfine;
 vector<double> slope;
 vector<double> exponent;
 vector<double> damp;
 vector<double> Q;
 vector<double> y_resp;
 vector<double> aa_psd;
 vector<double> vrs;
 vector<double> omega;
 vector<double> omega2;
 vector<double> omega4;
 vector<double> omegan;
 vector<double> omegan2;
 void initialize(void);
 void files(void);
 void frequency(void);
 void select_parameters(void);
 void header_read(long &number_lines,double &dt,double &sr);
 void rainflow(void);
 void line_test(int &kflag,const std::string & str);
 void transmitted_pressure_coefficients(double &fn,double &dam);
 void response(void);
 void SDOF_pressure_response(void);
 void output_files(void);
 void write_damage(void);
 void initialize_dr(void);
 void read_input_parameters(void);
 void generate_sample_psd(void);
 void generate_sample_psd2(void);
 void interp_sam(void);
 void fds_sample(void);
	 void compare(void);
	 void scalepsd(void);
	 void pressure_rms_sam(void);
	 void vrms_sam(void);
	 void drms_sam(void);
	 void writefile(void);
 void checklow(void);
 void low_core(void);
 void PRINTPSD(void);
 void writebest(void);
 void Dirlik_core(void);
 void interp(void);
 void cm0(void);
 void cm1(void);
 void cm2(void);
 void cm4(void);
 void weights_onethird(void);
 void convert_psd_to_spl(void);
};
int mflag; // leave outside class
void FDS::initialize(void)
{
 iq=0;
 np=0;
 numBytes=300;
 number_lines=0;
 pi=atan2(0.,-1.);
 tpi=2*pi;
 record =1.0e+90;
 pressure_rmslow=1.0e+50;
 vrmslow=1.0e+50;
 drmslow=1.0e+50;
 ratio_record=1.0e+50;
}
int main( int argc, char *argv[] )
{
 FDS fds;
 int mflag=0;
	if(argc> 1 )
	{
			printf("\n\n argc = %d ", argc );
			printf("\n argv[0] %s ", argv[0] );
			printf("\n argv[1] %s ", argv[1] );
			strcpy (fds.filename[0], argv[1] );
			mflag=1;
	}
	printf("\n envelope_fds_acoustic.cpp version 1.4, December 31, 2013 \n");
	printf("\n By Tom Irvine Email: tom@vibrationdata.com \n");
	printf("\n This program calculates an SPL to envelope an acoustic time ");
	printf("\n history in terms of fatigue damage spectra. The time history ");
	printf("\n may be nonstationary. \n");
 printf("\n Assume a uniform, fully-correlated acoustic field. \n");
 time_t start = time(0);
 fds.initialize();
	fds.files();
 fds.select_parameters();
	fds.header_read(fds.number_lines,fds.dt,fds.sr);
	fds.frequency();
 fds.read_input_parameters(); // Read input parameters from keyboard
//	printf(" out of frequency ");
	for(fds.iy=0;fds.iy<fds.damp.size();fds.iy++) { if(fds.iopt==1) { printf("\n\n calculating"); } fds.SDOF_pressure_response(); } printf("\n"); // printf("\n output files "); //************************************************************************** // printf("\n fds.ntrials %ld \n",fds.ntrials); double max=32767.; for(fds.ik=0; fds.ik<fds.ntrials; fds.ik++) { // printf(" %ld \n",fds.ik); if(rand()/max>0.5 &#124;&#124; fds.ik<20) { // printf("\n method 1 \n"); fds.generate_sample_psd(); // Generate the sample psd } else { // printf("\n method 2 \n"); fds.generate_sample_psd2(); } // printf(" ref 1 \n"); fds.interp_sam(); // Interpolate the sample psd at 1/24 octave // printf(" ref 2 \n"); fds.fds_sample(); // Calculate the fds of the sample psd // printf(" ref 3 \n"); fds.compare(); // Compare the sample fds with the reference fds // printf(" ref 4 \n"); fds.scalepsd(); // scale the psd // printf(" ref 5 \n"); fds.pressure_rms_sam(); // calculate the pressure_rms value // printf(" rms %8.4g \n",fds.pressure_rms); // printf(" ref 6 \n"); // fds.vrms_sam(); // printf(" ref 7 \n"); // fds.drms_sam(); // printf(" ref 8 \n"); // fds.checklow(); if( fds.ratio_trial_max < fds.ratio_record) { fds.low_core(); fds.ratio_record=fds.ratio_trial_max; } // printf(" ref 9 \n"); // printf("\n Trial: drms=%10.4g vrms=%10.4g pressure_rms=%10.4g ",fds.drms,fds.vrms,fds.pressure_rms); // printf("\n Record: drms=%10.4g vrms=%10.4g pressure_rms=%10.4g \n",fds.drmsp,fds.vrmsp,fds.pressure_rmsp); printf("\n Trial: ratio=%10.4g \t pressure_rms=%10.4g ",fds.ratio_trial_max,fds.pressure_rms); printf("\n Record: record=%10.4g \t pressure_rms_record=%10.4g \n",fds.ratio_record,fds.pressure_rmslow); // getch(); } printf("\n\n________________________________________________________________________"); for(long i=0; i < fds.nbreak; i++) { fds.f_sam[i]=fds.xf[i]; fds.apsd_sam[i]=fds.xapsd[i]; } fds.interp_sam(); // Interpolate the optimum psd at 1/24 octave fds.fds_sample(); // Calculate the fds of the optimum psd printf("\n\nOptimum Case\n"); fds.PRINTPSD(); fds.writebest(); // printf("\n drms=%10.4g vrms=%10.4g pressure_rms=%10.4g \n",fds.drmsp,fds.vrmsp,fds.pressure_rmsp); if(fds.iu==1) { printf("\n pressure_rms=%10.4g psi rms\n",fds.pressure_rmsp); printf("\n________________________________________________________________________"); printf("\n\n Optimum PSD written to file: best.psd - freq(Hz) pressure(psi^2/Hz) \n"); printf(" Optimum SPL written to file: best.spl - freq(Hz) SPL(dB) \n"); } else { printf("\n pressure_rms=%10.4g Pa rms\n",fds.pressure_rmsp); printf("\n________________________________________________________________________"); printf("\n\n Optimum PSD written to file: best.psd - freq(Hz) pressure(Pa^2/Hz) \n"); printf(" Optimum SPL written to file: best.spl - freq(Hz) SPL(dB) \n"); } printf("\n Overall SPL = %8.4g dB reference 20 micro Pa \n",fds.dB); fds.output_files(); fds.pFile[0]=fopen("vrs.txt","w"); for (long i=0;i<fds.aa_psd.size();i++) { fprintf(fds.pFile[0]," %12.5g \t %12.5g \n",fds.fn[i],fds.vrs[i]); } fclose(fds.pFile[0]); printf("\n The fatigue metric is: transmitted pressure (peak-valley)/2 \n"); if(fds.iu==1) { printf(" The reference area is 1 in^2 for comparison purposes only. \n"); } else { printf(" The reference area is 1 m^2 for comparison purposes only. \n"); } printf("\n\n Please call these files into your own graphics program. \n"); // printf("\n\n\n Input fds written to file: input.fds "); // printf("\n Optimum fds written to file: best.fds "); //************************************************************************** time_t end = time(0); double time = difftime(end, start); printf("\n Elapsed Time = %8.4g sec \n",time); printf("\n\n Press any key to exit.\n"); getch(); exit(1); } void FDS::output_files(void) { if(iopt==1){printf("\n");} printf("\n\n FDS Output files - fn(Hz) Damage Index \n\n"); for(i=0;i<iq;i++) { printf(" %s \n",input_fds_filename[i]); } printf("\n"); for(i=0;i<iq;i++) { printf(" %s \n",output_fds_filename[i]); } } void FDS::SDOF_pressure_response(void) { double ff; sprintf (temp_Q,"%g", Q[iy]); // printf(" Q=%g \n",Q[iy]); // printf(" %s \n",temp_Q); for (int count = 0; count < sizeof(temp_Q); count++) { if (temp_Q[count] == '.') { temp_Q[count] = 'p'; break; } } initialize_dr(); double dd; for(int ijk=0;ijk<fn.size();ijk++) { if(iopt==1) { printf("."); } // printf(" %ld %ld \n",ijk,fn.size()); ff=fn[ijk]; dd=damp[iy]; transmitted_pressure_coefficients(ff,dd); response(); // fprintf(pFile[1]," %14.7e \t %14.7e\n",fn[ijk],xabs[ijk]); } // fclose(pFile[1]); // printf("\n bb \n"); write_damage(); // printf("\n cc \n"); } void FDS::initialize_dr(void) { long long HEIGHT="0;" dr.resize(HEIGHT); for(i=0; i < HEIGHT; i++) { dr[i].resize(WIDTH); } std::vector<float> row;
 float a=0.;
 for(i=0;i<width;i++) { row.push_back(a); } for(i=0; i< HEIGHT;i++) { dr.push_back(row); } } void FDS::write_damage(void) { double ff; char temp_b[10]; for(i=0;i<exponent.size();i++) { sprintf (temp_b,"%g", exponent[i]); for (int count = 0; count < sizeof(temp_b); count++) { if (temp_b[count] == '.') { temp_b[count] = 'p'; break; } } //**** strcpy (filename[2],"input_damage_Q" ); strcat (filename[2],temp_Q); strcat (filename[2],"_b"); strcat (filename[2],temp_b ); strcat (filename[2],".txt" ); pFile[2] = fopen(filename[2], "w"); strcpy(input_fds_filename[iq],filename[2]); //**** strcpy (filename[4],"best_envelope_damage_Q" ); strcat (filename[4],temp_Q); strcat (filename[4],"_b"); strcat (filename[4],temp_b ); strcat (filename[4],".txt" ); strcpy(output_fds_filename[iq],filename[4]); //**** iq++; for(int ijk=0;ijk<fn.size();ijk++) { ff=fn[ijk]; fprintf(pFile[2]," %14.7e \t %14.7e\n",fn[ijk],dr[ijk][i]); if(iopt==2) { printf(" fn=%g Hz \t Q=%g \t b=%g \t Damage=%8.4g\n", fn[ijk],Q[iy],exponent[i],dr[ijk][i]); fprintf(pFile[3]," fn=%g Hz \t Q=%g \t b=%g \t Damage=%8.4g\n", fn[ijk],Q[iy],exponent[i],dr[ijk][i]); } reference_fds[iy][i][ijk]=dr[ijk][i]; // damp, exp, fn } fclose(pFile[2]); } } void FDS::frequency(void) { const long NUM=2000; oct=(1./24.); octh=oct/2; delta=pow(2,octh)-1./pow(2,octh); printf("\n t1=%9.5g sec dt=%8.4g sec\n",t1,dt); printf("\n\n number of input points = %ld \n\n",np); if(dt<1.0e-10){exit(1);} xabs.push_back(0.); double fr; j=1; while(1) { xabs.push_back(0.); fr= fn[0]*pow(2.,j*oct); if( fr> sr/2. )
		{
 fn.push_back(fr);
		 break;
 }
 else
 {
 if( fr>=flast)
 {
 fn.push_back(flast);
 break;
 }
 else
 {
 fn.push_back(fr);
 }
 }
 j++;
 }
 jnum=j;
 n_ref=j;
 double ofr;
 for(long j=0;j<fn.size();j++) { ofr=tpi*fn[j]; omega.push_back(ofr); omegan.push_back(ofr); omega2.push_back(pow(ofr,2)); omegan2.push_back(pow(ofr,2)); omega4.push_back(pow(ofr,4)); } /* for(long j=0;j<fn.size();j++) { printf(" j=%ld fn=%8.4g flast=%8.4g \n",j,fn[j],flast); } printf("\n fn.back %8.4g \n",fn.back()); getch(); */ } void FDS::select_parameters(void) { double ff; /********************* intialize frequency array *********************/ // printf("\n Select rainflow dimension: "); // printf("\n 1=acceleration 2=relative displacement \n"); // scanf("%d",&iresp); iresp=1; // printf("\n Select natural frequency analysis options: "); // printf("\n 1=spectrum 2=single fn \n"); // scanf("%d",&iopt); iopt=1; if(iopt==1) { printf("\n\n Enter starting frequency (Hz) \n"); scanf("%lf",&ff); fn.push_back(ff); printf("\n\n Enter ending frequency (Hz) \n"); scanf("%lf",&flast); } pFile[3]=fopen(filename[3],"w"); printf("\n Enter the number of amplification factor Q cases (max=3): \n"); scanf("%d",&iqcases); if(iqcases>3){iqcases=3;}
 double qq;
	for(i=0;i<iqcases;i++) { printf("\n Enter amplification factor %ld: ",i+1); scanf("%lf",&qq); Q.push_back(qq); damp.push_back(1./(2.*Q[i])); } printf("\n Enter the number of fatigue exponent cases (max=3): \n"); scanf("%d",&fcases); if(fcases>3){fcases=3;}
 double ee;
 for (i=0;i<fcases;i++) { printf("\n Enter fatigue exponent: %ld \n",i+1); scanf("%lf",&ee); exponent.push_back(ee); } } void FDS::response(void) { /********************** response calculation *******************/ // printf(" %d fn=%8.4g Hz \n ",index,fn[index]); double yy; double xb=0.; double xbb=0.; double yb=0.; double ybb=0.; for(long i=0;i<np;i++) { x=0.; yy=yin[i]; x= a1*xb +a2*xbb +af1*yy +af2*yb +af3*ybb; if(fabs(x)> xabs[index]){xabs[index]=fabs(x);}
	 xbb=xb;
	 xb= x;
 ybb=yb;
	 yb=yy;
 y_resp.push_back(x);
 }
// printf(" xabs[index]=%8.4g ",xabs[index]);
// printf(" fn=%8.4g \n",fn[index]);
 rainflow();
 y_resp.clear();
}
void FDS::rainflow(void)
{
//**************************************************************
// printf("\n NP=%ld \n",NP);
 long nrows=1200000;
 DYNAMIC_MATRIX(B,nrows,2);
 float C[100];
 float AverageMean[100];
 float MaxPeak[100];
 float MinValley[100];
 float MaxAmp[100];
 float AverageAmp[100];
 long k;
 long last_a=0;
 vector<float> a;
 double slope1;
 double slope2;
 double X;
 double Y;
 double ymax;
 double amax=0;
 for(long i=0;i<y_resp.size();i++) { if(fabs(y_resp[i])>amax)
	 {
	 amax=fabs(y_resp[i]);
	 }
	}
// printf("\n amax=%8.4g \n",amax);
	k=0;
//	a[k]=y[k];
 a.push_back(y_resp[k]);
	k=1;
 long NP=y_resp.size();
 for(long i=1;i<(np-1);i++) { slope1=( y_resp[i]-y_resp[i-1]); slope2=(y_resp[i+1]-y_resp[i]); if((slope1*slope2)<=0) { a.push_back(y_resp[i]); k++; } } long hold=last_a; a.push_back(y_resp.back()); k++; last_a=k-1; long mina=100000; long maxa=-mina; for(long i=0;i<=last_a;i++) { if(a[i]<mina) { mina=a[i]; } if(a[i]>maxa)
		{
			maxa=a[i];
		}
//		printf(" %8.4g \n",a[i]);
	}
	long num=long(maxa-mina)+1;
	long n=0;
	long i=0;
	long j=1;
	long ijk=0;
	float sum=0;
	long kv=0;
	long LLL=last_a;
 long nkv=0;
// printf("\n percent completed \n");
	while(1)
	{
 	Y=(fabs(a[i]-a[i+1]));
 	X=(fabs(a[j]-a[j+1]));
 	if(X>=Y && Y>0)
		{
 if(Y>ymax)
 {ymax=Y;}
 	if(i==0)
			{
 		n=0;
 		sum+=0.5;
// B[kv][3]=a[i+1];
// B[kv][2]=a[i];
 B[kv][1]=0.5;
 B[kv][0]=Y;
 kv++;
/*
				for(k=0;k<last_a;k++) { a[k]=a[k+1]; } */ a.erase (a.begin()); last_a--; i=0; j=1; } else { sum+=1; // B[kv][3]=a[i+1]; // B[kv][2]=a[i]; B[kv][1]=1.; B[kv][0]=Y; kv++; n=0; a.erase (a.begin()+(i+1)); a.erase (a.begin()+i); last_a-=2; i=0; j=1; } nkv++; if(nkv==3000) { double ratio = fabs((last_a)/double(LLL)); // printf(" %3.1lf \n",(1-ratio)*100.); nkv=0; } } else { i++; j++; } if((j+1)>last_a)
		{
 	break;
		}
	}
	for(i=0;i<(last_a);i++) { Y=(fabs(a[i]-a[i+1])); if(Y>0)
		{
 	sum+=0.5;
 // B[kv][3]=a[i+1];
 // B[kv][2]=a[i];
 B[kv][1]=0.5;
 B[kv][0]=Y;
 kv++;
		}
 if(Y>ymax)
 {ymax=Y;}
	}
 std::vector<float> row;
 float nc;
 for(j=0;j<exponent.size();j++) { damage=0; for(i=0;i<=kv;i++) { Y=B[i][0]/2.; // amp nc=B[i][1]; // cycle damage+=nc*pow(Y,exponent[j]); } row.push_back(damage); } dr.push_back(row); DELETE_MATRIX(B,nrows); } void FDS::files(void) { printf("\n Select pressure unit: \n 1=psi 2=Pa \n"); scanf("%d",&iu); /******************************** open files *************************/ if(mflag==0) { printf( "\n\n The applied pressure file must contain two columns: \n"); if(iu==1) { printf( " time(sec) pressure(psi) \n"); } else { printf( " time(sec) pressure(Pa) \n"); } printf( "\n 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"); } } void FDS::header_read(long &number_lines,double &dt,double &sr) { double ymax=0.; long numBytes=300; char str[numBytes]; long j=0; int kflag; double tb=-1.0e+90; double dtmin=1.0e+90; double dtmax=-dtmin; float ya; printf("\n"); np=0.; while(fgets(string_a,numBytes,pFile[0] )>0)
	{
 line_test(kflag,string_a);
		if(kflag==1)
		{
			printf("H: %s ",string_a);
		}
		else
		{
			i=number_lines;
			if(strrchr(string_a,',' ) != NULL )
			{
				sscanf(string_a,"%lf, %f", &t,&ya );
				yin.push_back(ya*1); // convert pressure to pressure
 j++;
 np++;
			}
			else
			{
				sscanf(string_a,"%lf %f", &t,&ya ); // convert pressure to pressure
 yin.push_back(ya*1);
 if(yin[j]>ymax)
 {
 ymax=yin[j];
 }
 j++;
 np++;
			}
 if(number_lines==0)
			{
				t1=t;
			}
			if(t<tb && number_lines>=1)
			{
				printf("\n Error: time values must ascend. \n");
				printf("\n %s \n",string_a);
			 printf("\n t=%11.7e \t tb=%11.7e \n",t,tb);
				printf("\n Press any key to exit \n");
				getch();
				exit(1);
			}
			if(number_lines>=1)
			{
				if(dtmin>(t-tb))
				{
					dtmin=(t-tb);
				}
				if(dtmax<(t-tb)) { dtmax=(t-tb); } } tb=t; number_lines++; if(number_lines==MAX) { printf("\n\n ** warning: maximum line number reached. ** \n"); break; } } } rewind(pFile[0]); dt=(t-t1)/double(number_lines-1); sr=1./dt; printf("\n Time span: %10.5g to %10.5g sec \n",t1,t); printf("\n dtmin= %11.6g sec ",dtmin); printf("\n dt = %11.6g sec ",dt); printf("\n dtmax= %11.6g sec \n",dtmax); printf("\n srmin= %11.6g samples/sec ",1./dtmax); printf("\n sr = %11.6g samples/sec ",sr); printf("\n srmax= %11.6g samples/sec \n",1./dtmin); printf("\n ymax=%8.4g \n",ymax); if( dtmin < 0.99 *dtmax) { int nsr; printf("\n *** Time Step Warning *** \n"); printf("\n dtmax = %8.4g dtmin = %8.4g \n",dtmax,dtmin); printf("\n The time step must be constant. \n"); printf("\n Enter new sample rate? 1=yes 2=no \n"); scanf("%d",&nsr); if(nsr==1) { printf("\n\n Assume constant sample rate. Enter sample rate (samples/sec) \n"); scanf("%lf",&sr); dt=1./sr; printf("\n\n sr = %11.6g samples/sec ",sr); printf("\n Nyquist Freq = %lf Hz \n",sr/2.); printf("\n dt = %11.6g sec \n",dt); } else { printf("\n Press any key to exit \n"); getch(); exit(1); } } } void FDS::line_test(int &kflag,const std::string & str) { // detect header lines // printf("j=%ld. %s\n",j,string); kflag=0; 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 FDS::transmitted_pressure_coefficients(double &fn,double &dam)
{
 double C,E,K,S,Sp;
 double omega=2.*pi*fn;
 double omegad=omega*sqrt(1.-pow(dam,2));
 E=exp(-dam*omega*dt);
 K=omegad*dt;
 C=E*cos(K);
 S=E*sin(K);
 Sp=S/K;
 a1=2*C;
 a2=-pow(E,2);
 af1=1.-Sp;
 af2=2.*(Sp-C);
 af3=pow(E,2)-Sp;
}
void FDS::read_input_parameters(void)
{
/*
	printf("\n\n Enter goal. Minimize: \n");
	printf("\n 1=acceleration ");
	printf("\n 2=acceleration, velocity ");
	printf("\n 3=acceleration, velocty, displacement ");
	printf("\n 4=acceleration, displacement ");
	printf("\n 5=displacement \n");
 scanf("%ld", &goal);
*/
 goal=1;
	printf("\n\n Input number of trials ");
 scanf("%ld", &ntrials);
/*
	printf("\n Input number of trials for phase 2 (refinement phase): ");
 scanf("%ld", &ntrials2);
*/
 long MAX2=100000;
 if(ntrials> MAX2)
	{
		printf("\n The maximum number of trials is 100,000 \n");
		ntrials=MAX2;
	}
 printf("\n\n Enter number of breakpoints for envelope (min=3 max=12): ");
 scanf("%ld",&nbreak);
 if(nbreak < 3 )nbreak=3; if(nbreak> 12 )nbreak=12;
 printf("\n\n Constrain slopes? 1=yes 2=no ");
 scanf("%ld",&ic);
 if(ic==1)
 {
	 printf("\n Enter absolute slope limit (dB/octave): ");
	 scanf("%lf",&slopec);
 }
 else
 {
	 slopec = 60.;
 }
 slopec=fabs(slopec);
 slopec=(slopec/10.)/log10(2.);
// printf("\n\n n limit = %lf \n",slopec);
 printf("\n\n Further constrain initial slope to be positive ? 1=yes 2=no ");
 scanf("%ld",&initial);
 printf("\n\n Further constrain final slope to be negative ? 1=yes 2=no ");
 scanf("%ld",&final);
 printf("\n Enter the PSD duration (sec) ");
 scanf("%lf",&duration);
}
void FDS::generate_sample_psd()
{
// printf(" start generate sample psd \n");
 double max=32767.;
 f_sam.clear();
	long ijk;
	double temp;
//	 printf(" nbreak=%ld \n",nbreak);
	 // printf(" n_ref=%ld \n",n_ref);
	f_sam.push_back(fn[0]);
	f_sam.push_back(flast);
// printf(" continue generate sample psd \n");
 // generate some random numbers for frequency
 int iflag;
 int ia[20];
 ia[0]=0;
 ia[1]=fn.size()-1;
	 for(ijk=0; ijk < 10000; ijk++) { // printf(" ijk=%ld \n",ijk); for(i=2; i<nbreak; i++) { while(1) { iflag=0; index = long( n_ref*rand()/max ); // printf(" index=%ld \n",index); for(int nv=0;nv<i;nv++) { // printf(" ia[%d]=%d \n",nv,ia[nv]); if(index==ia[nv]) { iflag=1; break; } } if(iflag==0) { ia[i]=index; break; } // printf(" iflag=%d \n",iflag); } f_sam.push_back(fn[index]); } // sort the frequencies FSORT } // printf("\n **** %8.4g %8.4g **** \n",f_sam.back(),flast); // printf(" ref b \n"); // generate some random numbers for amplitude double apsd_max=0.; for(i=0; i< nbreak; i++) { // printf(" f_sam[%ld]=%8.4g \n",i,f_sam[i]); apsd_sam[i]=rand()/max; if(apsd_sam[i]>apsd_max)
		 {
		 apsd_max=apsd_sam[i];
		 }
	 }
	 for(i=0; i< nbreak; i++) { apsd_sam[i]*=0.001/apsd_max; } // printf(" ref c \n"); double amin=1.0e-12; for(i=0; i< nbreak; i++) { if(apsd_sam[i] < amin){apsd_sam[i]=amin;} if(apsd_sam[i]> (1./amin)){apsd_sam[i]=(1./amin);}
	 }
// printf(" ref d \n");
 // absolute slopes are limited to slopec
	ABS_SLOPE
	apsd_sam_last=apsd_sam[nbreak-1];
	printf("\n\n Phase 1, Trial %ld, PSD Coordinates \n",ik);
/*
 for(i=0; i< nbreak; i++) { printf(" f_sam=%8.4g apsd_sam=%8.4g \n",f_sam[i],apsd_sam[i]); } */ // printf("\n press any key \n"); // getch(); } void FDS::generate_sample_psd2() { f_sam.clear(); double max=32767.; double temp; f_sam.push_back(fn[0]); f_sam.push_back(flast); double aaa,bbb; if( rand()/max>0.2)
	{
		bbb=pow((rand()/max),3.);
	}
	else
	{
		bbb=0.1*rand()/max;
	}
	aaa=1.-bbb/2.;
	for(i=1; i<= (nbreak-2); i++) { f_sam.push_back(xf[i]*(aaa+bbb*rand()/max)); } // sort the frequencies FSORT FSPACE f_sam.back()=flast; // printf("\n **** %8.4g %8.4g **** \n",f_sam.back(),fn.back()); // generate some random number for amplitude // double db; bbb=pow((rand()/max),3.); aaa=1.-bbb/2.; for(i=0; i< nbreak; i++) { apsd_sam[i]=xapsd[i]*(aaa+bbb*rand()/max); if(apsd_sam[i]>2000.)
			{
				printf("\n Error: %ld %8.4g %8.4g \n",i,apsd_sam[i],xapsd[i]);
				exit(1);
			}
		}
 // add constraint
	ABS_SLOPE
 printf("\n\n Phase 2, Trial %ld, PSD Coordinates \n",ik);
//	PRINTPSD
}
void FDS::interp_sam(void)
{
 ppsd_samfine.clear();
	ppsd_samfine.push_back(apsd_sam[0]);
	for(i=0; i< (nbreak-1);i++) { slope.push_back(log(apsd_sam[i+1]/apsd_sam[i])/log(f_sam[i+1]/f_sam[i])); } for(i=1; i < fn.size(); i++) { for(j=0; j < (nbreak-1); j++) { if( ( fn[i]>= f_sam[j] ) && ( fn[i] <= f_sam[j+1] ) ) { ppsd_samfine.push_back(apsd_sam[j]*pow( ( fn[i] / f_sam[j] ), slope[j] )); break; } } } /* FILE *ppFile[1]; char pfilename[1][FILENAME_MAX]; ppFile[1] = fopen("psd.inp", "w"); for(i=0; i < ppsd_samfine.size(); i++) { fprintf(ppFile[1]," %8.4g \t %8.4g \n",fn[i],ppsd_samfine[i]); printf(" %8.4g \t %8.4g \n",fn[i],ppsd_samfine[i]); } fclose(ppFile[1]); */ // printf("\n press any key \n"); // getch(); } void FDS::fds_sample(void) { // printf("\n inside fds_sample \n"); double suma; double fn2; double d; double df; double tdamp; double tdamp2; double rho; double tdr; double c1,c2; double trans; // fn // ppsd_samfine // printf("\n"); // printf(" damp.size()=%ld \n",damp.size()); // printf(" fn.size()=%ld \n",fn.size()); // printf(" fn.size()=%ld \n",fn.size()); // printf(" exponent.size()=%ld \n",exponent.size()); for(int ijk=0;ijk<damp.size();ijk++) { // printf(" ijk=%ld \n",ijk); tdamp=2*damp[ijk]; tdamp2=pow(tdamp,2.); vrs.clear(); for(long i=0;i<fn.size();i++) // keep long { // printf(" n_ref=%ld \n",n_ref); aa_psd.clear(); for(long j=0;j<fn.size();j++) // keep long { rho = fn[j]/fn[i]; tdr=tdamp*rho; c1= pow(tdr,2.); c2= pow((1.- pow(rho,2.)),2.); trans=(1.+ c1 ) / ( c2 + c1 ); aa_psd.push_back(ppsd_samfine[j]*trans); // printf(" fn[j]=%8.4g c1=%8.4g c2=%8.4g \n",fn[j],c1,c2); // printf(" aa_psd=%8.4g ppsd_samfine[%ld]=%8.4g trans=%8.4g\n",aa_psd,j,ppsd_samfine[j],trans); // printf(" fn=%8.4g f_s=%8.4g aa_psd=%8.4g df=%8.4g\n",fn[i],fn[j],aa_psd,df); } //*************************************************************************************** cm0(); cm1(); cm2(); cm4(); //*************************************************************************************** EP=sqrt(m4/m2); ft_rms=sqrt(m0); vrs.push_back(ft_rms); /* if(i==15) { printf("\n fn=%8.4g m0=%8.4g m1=%8.4g pressure_rms=%8.4g \n",fn[i],m0,m1,pressure_rms); } */ for(long j=0;j<exponent.size();j++) // keep long { // printf(" j=%ld go to Dirlik \n",j); Dirlik_core(); // printf(" out of Dirlik \n"); d=0.; for(long jk=0;jk<amp.size();jk++) { d+=pow(amp[jk],exponent[j]); // printf(" d=%8.4g amp=%8.4g exponent=%8.4g \n",d,amp[jk],exponent[j]); } // printf(" ijk=%ld j=%ld i=%ld d=%8.4g amp.size=%ld \n",ijk,j,i,d,amp.size()); fatigue_damage[ijk][j][i]=d; } // printf("ref aaa press any key\n"); // getch(); } // fn } // damp } void FDS::Dirlik_core() { // 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*ft_rms; // printf(" rms %8.4g ",pressure_rms); 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; 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()); long num=long(0.5+cumu.back()); // printf("\n num %ld \n",num); xq.clear(); for(long i=0;i<num;i++) { xq.push_back(i); } // 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 \n",amp.size(),maxamp);
// printf(" end DDD");
}
void FDS::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<xq.size();i++) // keep long { peak_range.push_back(0.); } for(long i=1; i<xq.size(); i++) { ttt = xq[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 FDS::compare(void) { // printf("\n in compare\n"); double sc; scale=0; // printf(" %ld \n",damp.size()); // printf(" %ld \n",exponent.size()); // printf(" %ld \n",fn.size()); for(int ijk=0;ijk<damp.size();ijk++) { for(long j=0;j<exponent.size();j++) // keep long { for(long i=0;i<fn.size();i++) { // printf(" rfds=%8.4g fds=%8.4g \n",reference_fds[ijk][j][i],fatigue_damage[ijk][j][i]); sc=reference_fds[ijk][j][i]/fatigue_damage[ijk][j][i]; sc=pow(sc,(1./exponent[j])); if(sc>scale)
 {
 scale=sc;
 }
 }
 }
 }
// printf(" scale=%8.4g \n",scale);
}
void FDS::scalepsd(void)
{
	scale=pow(scale,2.);
// printf("\n scale = %8.4g nbreak=%ld \n",scale,nbreak);
	for(i=0; i< nbreak; i++) { // printf(" apsd_sam = %8.4g \n",apsd_sam[i]); apsd_sam[i]*=scale; } apsd_sam_last=apsd_sam[nbreak-1]; PRINTPSD(); //********************************************************************************* double fdd; double ratio; double sce; ratio_trial_max=0.; for(int ijk=0;ijk<damp.size();ijk++) { for(long j=0;j<exponent.size();j++) // keep long { sce=pow(scale,exponent[j]/2.); for(long i=0;i<fn.size();i++) { // printf(" rfds=%8.4g fds=%8.4g \n",reference_fds[ijk][j][i],fatigue_damage[ijk][j][i]); fdd=fatigue_damage[ijk][j][i]*sce; ratio=fdd/reference_fds[ijk][j][i]; if(ratio>ratio_trial_max)
 {
 ratio_trial_max=ratio;
 }
 }
 }
 }
//*********************************************************************************
}
void FDS::cm0(void)
{
// m0=m0+aa_psd*df;
 double ra=0.;
 for(long i=0; i < (aa_psd.size()-1); i++) { s=( log( aa_psd[i+1]/aa_psd[i] )/log( fn[i+1]/fn[i] ) ); if(s < -1.0001 &#124;&#124; s> -0.9999 )
 {
 ra+= ( aa_psd[i+1] * fn[i+1]- aa_psd[i]*fn[i])/( s+1.);
 }
 else
 {
 ra+= aa_psd[i]*fn[i]*log( fn[i+1]/fn[i]);
 }
 }
 m0=ra;
}
void FDS::cm1(void)
{
// m1=m1+aa_psd*fn[j]*df;
 double ra=0.;
 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++) { 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 FDS::cm2(void)
{
// m2=m2+aa_psd*pow(fn[j],2)*df;
 double ra=0.;
 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++) { 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 FDS::cm4(void)
{
// m4=m4+aa_psd*pow(fn[j],4)*df;
 double ra=0.;
 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++) { 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 FDS::pressure_rms_sam(void)
{
 double ra=0.;
 vector<double> s;
 pressure_rms=0.;
	 for(long i=0; i < nbreak-1; i++) { s.push_back( log( apsd_sam[i+1]/apsd_sam[i] )/log( f_sam[i+1]/f_sam[i] ) ); if(s[i] < -1.0001 &#124;&#124; s[i]> -0.9999 )
			 {
 ra+= ( apsd_sam[i+1] * f_sam[i+1]- apsd_sam[i]*f_sam[i])/( s[i]+1.);
			 }
		 else
			 {
 ra+= apsd_sam[i]*f_sam[i]*log( f_sam[i+1]/f_sam[i]);
			 }
		 }
 pressure_rms=sqrt(ra);
}
void FDS::writefile(void)
{
	for( i=0; i < nbreak; i++) { fprintf(pFile[2],"%12.5e %12.5e\n",xf[i],xapsd[i]); f_sam[i]=xf[i]; apsd_sam[i]=xapsd[i]; } /* for( i=0; i < n_ref; i++) { fprintf(pFile[3],"%12.5e %12.5e\n",xffine[i],xfds[i]); } */ } void FDS::PRINTPSD(void) { double s1,s2; printf("\n Freq PSD Slope Slope \n (Hz) (G^2/Hz) (N) (dB/octave) \n"); for(i=0; i<nbreak; i++) { if( i < (nbreak-1)) { s1=log(apsd_sam[i+1]/apsd_sam[i])/log(f_sam[i+1]/f_sam[i]); s2=s1*10.*log10(2.); printf("%12.2f %12.5f %12.4f %12.4f\n",f_sam[i],apsd_sam[i],s1,s2 ); } else { printf("%12.2f %12.5f\n",f_sam[i],apsd_sam[i],s1,s2 ); } } } void FDS::writebest() { double ref; if(iu==1) { ref=2.9e-09; } else { ref=20.e-06; } strcpy(filename[2],"best.psd"); pFile[2]=fopen(filename[2],"w"); strcpy(filename[5],"best.spl"); pFile[5]=fopen(filename[5],"w"); for(long i=0; i < nbreak; i++) { fprintf(pFile[2],"%12.5e \t %12.5e\n",xf[i],xapsd[i]); f_sam[i]=xf[i]; apsd_sam[i]=xapsd[i]; } fclose(pFile[2]); convert_psd_to_spl(); double rms=0.; double spm2; for(long i=0; i <= fc_num; i++) { spl[i]=0.; // printf(" i=%ld sum_pressure=%8.4g \n",i,sum_pressure[i]); if(sum_pressure[i]>1.0e-20)
 {
 spm2=sqrt(sum_pressure[i]);
 spl[i]=20.*log10(spm2/ref);
 fprintf(pFile[5],"%12.5e \t %12.5e\n",fc[i],spl[i]);
// printf(" i=%ld spm2=%8.4g spl=%8.4g ref=%8.4g \n",i,spm2,spl[i],ref);
 rms+=sum_pressure[i];
 }
 }
 fclose(pFile[5]);
 dB=20.*log10(sqrt(rms)/ref);
 int kj=0;
 for(int ijk=0;ijk<damp.size();ijk++) { for(long j=0;j<exponent.size();j++) // keep long { pFile[4]=fopen(output_fds_filename[kj], "w"); kj++; for(long i=0;i<fn.size();i++) { fprintf(pFile[4],"%12.5e \t %12.5e\n",fn[i],fatigue_damage[ijk][j][i]); } fclose(pFile[4]); } } } void FDS::low_core() { // printf("\n in low core \n"); pressure_rmslow=pressure_rms; pressure_rmsp=pressure_rms; ikbest=ik; nnn=ntrials; for( i=0; i < nbreak; i++) { xf[i]=f_sam[i]; xapsd[i]=apsd_sam[i]; // printf(" xf=%8.4g \t xapsd=%8.4g \n",xf[i],xapsd[i]); } } void FDS::weights_onethird(void) { octbw=( pow(2.,(1./6.)) - 1./(pow(2.,(1./6.))) ); fc[0]=10.; fc[1]=12.5; fc[2]=16.; fc[3]=20.; fc[4]=25.; fc[5]=31.5; fc[6]=40.; fc[7]=50.; fc[8]=63.; fc[9]=80.; fc[10]=100.; fc[11]=125.; fc[12]=160.; fc[13]=200.; fc[14]=250.; fc[15]=315.; fc[16]=400.; fc[17]=500.; fc[18]=630.; fc[19]=800.; fc[20]=1000.; fc[21]=1250.; fc[22]=1600.; fc[23]=2000.; fc[24]=2500.; fc[25]=3150.; fc[26]=4000.; fc[27]=5000.; fc[28]=6300.; fc[29]=8000.; fc[30]=10000.; fc[31]=12500.; fc[30]=16000.; fc[31]=20000.; fc[32]=25000.; fc[33]=31500.; fc_num=33; } void FDS::convert_psd_to_spl(void) { double df; delta=pow(2,octh)-1./pow(2,octh); weights_onethird(); for(i=0; i<=fc_num; i++) { sum_pressure[i]=0.; fband[i]=fc[i]/(pow(2.,(1./6.))); } for(i=0; i<fn.size(); i++) { df=fn[i]*delta; ppsd_samfine[i]*=df; for(j=0; j<fc_num; j++) { if( fn[i]>= fband[j] && fn[i] < fband[j+1] )
			{
			 sum_pressure[j]+=1*ppsd_samfine[i];
			 break;
			}
 }
	}
}
</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/matlab2/envelope_fds_acoustic.cpp">(-&gt;オリジナル)</a>
/ <label>アドレス: <input type="text" name="naked_post_url" value="http://www.vibrationdata.com/matlab2/envelope_fds_acoustic.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>