/*========================*/ /* gccEQSP.c */ /*========================*/ #include #include #include #include #define Nmax 32768 #define Mmax 200 void main_part(char fnameR[],char fnameW[]); void ERES(int nn,double dt,double h,double ddy[Nmax],int nt, double tk[Mmax],double resacc[Mmax],double resvel[Mmax],double resdis[Mmax]); void IACC(int nn,double dt,double ddy[Nmax],double dy[Nmax],double y[Nmax], double *ddymax,double *dymax,double *ymax); void CRAC(int nn,double dt,double ddy[Nmax],double dy[Nmax],double y[Nmax], double ddymax,double dymax,double ymax); /*---------------------------------------------------------------------------*/ int main(int argc,char *argv[]) { char fnameR[50]; /*入力ファイル名*/ char fnameW[50]; /*出力ファイル名*/ strcpy(fnameR,argv[1]); strcpy(fnameW,argv[2]); main_part(fnameR,fnameW); return 0; } /*---------------------------------------------------------------------------*/ void main_part(char fnameR[],char fnameW[]) { FILE *fin,*fout; char dat[256]; int i; char scom[256]; /*コメント*/ double dt; /*サンプリングピッチ(sec)*/ int nn; /*サンプリングデータ数*/ double ddy[Nmax]; /*加速度時刻歴*/ double dy[Nmax]; /*速度時刻歴*/ double y[Nmax]; /*変位時刻歴*/ double h=0.05; /*減衰定数*/ int nt=Mmax; /*周期計算個数*/ double tk[Mmax]; /*計算周期*/ double resacc[Mmax]; /*加速度応答*/ double resvel[Mmax]; /*速度応答*/ double resdis[Mmax]; /*変位応答*/ double ddymax,dymax,ymax; /*データ入力*/ fin=fopen(fnameR,"r"); fgets(dat,sizeof dat,fin);strcpy(scom,dat); fgets(dat,sizeof dat,fin);strtok(dat,",");dt=atof(strtok(NULL,",")); fgets(dat,sizeof dat,fin);strtok(dat,",");nn=atoi(strtok(NULL,",")); for(i=0;i<=nn-1;i++){ fgets(dat,sizeof dat,fin);ddy[i-1]=atof(strtok(dat,",")); } fclose(fin); /*計算周期(sec)設定*/ for(i=0;i<=nt-1;i++){ tk[i]=pow(10.0,log10(2.0*dt)+(log10(10.0)-log10(2.0*dt))/(double)nt*(double)i); } /*応答スペクトル計算*/ IACC(nn,dt,ddy,dy,y,&ddymax,&dymax,&ymax); /*加速度時刻歴数値積分*/ CRAC(nn,dt,ddy,dy,y,ddymax,dymax,ymax); /*加速度時刻歴基線補正*/ ERES(nn,dt,h,ddy,nt,tk,resacc,resvel,resdis); /*応答スペクトル計算*/ /*データ書き込み*/ fout=fopen(fnameW,"w"); fprintf(fout,"%s",scom); fprintf(fout,"nt,%d\n",nt); fprintf(fout,"period(sec),acc(gal),vel(kine),dis(cm)\n"); for(i=0;i<=nt-1;i++){ fprintf(fout,"%e,%e,%e,%e\n",tk[i],resacc[i],resvel[i],resdis[i]); } fclose(fout); } /*---------------------------------------------------------------------------*/ void ERES(int nn,double dt,double h,double ddy[Nmax],int nt, double tk[Mmax],double resacc[Mmax],double resvel[Mmax],double resdis[Mmax]) { /*応答値算出 */ /* nn :加速度時刻歴総数 */ /* dt :時間間隔(入力値) */ /* h :減衰定数(入力値) */ /* ddy[] :加速度時刻歴(入力値) */ /* nt :計算周期総数[データ格納0〜nt-1](入力値)*/ /* tk[] :計算周期(入力値) */ /* resacc[] :加速度応答スペクトル(計算値) */ /* resvel[] :速度応答スペクトル(計算値) */ /* resdis[] :変位応答スペクトル(計算値) */ int i,m; double accmax,velmax,dismax; /*応答最大値(gal,kine,cm)*/ double w,w2,hw,wd,wdt,e,cwdt,swdt; double ss,cc,s1,c1,s2,c2,s3,c3; double A11,A12,A21,A22,B11,B12,B21,B22; double dxf,xf,ddym,ddyf,x,dx,ddx; for(i=0;i<=nt-1;i++){ w=2.0*M_PI/tk[i]; w2=w*w; hw=h*w; wd=w*sqrt(1.0-h*h); wdt=wd*dt; e=exp(-hw*dt); cwdt=cos(wdt); swdt=sin(wdt); A11=e*(cwdt+hw*swdt/wd); A12=e*swdt/wd; A21=-e*w2*swdt/wd; A22=e*(cwdt-hw*swdt/wd); ss=-hw*swdt-wd*cwdt; cc=-hw*cwdt+wd*swdt; s1=(e*ss+wd)/w2; c1=(e*cc+hw)/w2; s2=(e*dt*ss+hw*s1+wd*c1)/w2; c2=(e*dt*cc+hw*c1-wd*s1)/w2; s3=dt*s1-s2; c3=dt*c1-c2; B11=-s2/wdt; B12=-s3/wdt; B21=(hw*s2-wd*c2)/wdt; B22=(hw*s3-wd*c3)/wdt; accmax=2.0*hw*fabs(ddy[0])*dt; velmax=fabs(ddy[0])*dt; dismax=0.0; dxf=-ddy[0]*dt; xf=0.0; for(m=1;m<=nn-1;m++){ ddym=ddy[m]; ddyf=ddy[m-1]; x=A12*dxf+A11*xf+B12*ddym+B11*ddyf; dx=A22*dxf+A21*xf+B22*ddym+B21*ddyf; ddx=-2.0*hw*dx-w2*x; if(accmax<=fabs(ddx))accmax=fabs(ddx); if(velmax<=fabs(dx))velmax=fabs(dx); if(dismax<=fabs(x))dismax=fabs(x); dxf=dx; xf=x; } resacc[i]=accmax; resvel[i]=velmax; resdis[i]=dismax; } } /*---------------------------------------------------------------------------*/ void IACC(int nn,double dt,double ddy[Nmax],double dy[Nmax],double y[Nmax], double *ddymax,double *dymax,double *ymax) { /*加速度時刻歴数値積分 */ /* nn :加速度時刻歴データ総数 */ /* dt :時間間隔(入力値) */ /* ddy[] :加速度時刻歴(入力値) */ /* dy[] :速度時刻歴(計算出力) */ /* y[] :変位時刻歴(計算出力値) */ /* ddymax :加速度最大値(計算出力値)*/ /* dymax :速度最大値(計算出力値) */ /* ymax :変位最大値(計算出力値) */ int i; *ddymax=0.0; *dymax=0.0; *ymax=0.0; dy[0]=0.0; y[0]=0.0; for(i=1;i<=nn-1;i++){ dy[i]=dy[i-1]+(ddy[i-1]+ddy[i])*dt/2.0; y[i]=y[i-1]+dy[i-1]*dt+(ddy[i-1]/3.0+ddy[i]/6.0)*dt*dt; if(*ddymax