%該程序為vidic主程序
clear;%清除工作空間變量；
clc;format compact;format long g;
%per=[0:0.2:1.0 1.5:0.5:2 3:1:6];%設定自振周期范圍,分別采用不同的密度；
per=0:0.2:6;
sp=length(per);
%采用不同的屈服強度系數；
% ebx(1)=0.167;
% ebx(2)=0.2;
% ebx(3)=0.25;
ebx(4)=0.333;%%Ru=3,強度折減系數；
% ebx(5)=0.5;
for iy=1:1
	%=加載地震波=-------------------------------------------------------------------
	%next one 46
	fida(1)=fopen('d:/berkleywave/northridge/nor_chl160_a.txt','r+');
	fidv(1)=fopen('d:/berkleywave/northridge/nor_chl160_v.txt','r+');
	ddt(1)=0.01;
	%next one 108
	fida(2)=fopen('d:/berkleywave/northridge/nor_jen292_a.txt','r+');
	fidv(2)=fopen('d:/berkleywave/northridge/nor_jen292_v.txt','r+');
	ddt(2)=0.02;
	%next one 124
	fida(3)=fopen('d:/berkleywave/northridge/nor_ldm334_a.txt','r+');
	fidv(3)=fopen('d:/berkleywave/northridge/nor_ldm334_v.txt','r+');
	ddt(3)=0.005;
	%next one 175
	fida(4)=fopen('d:/berkleywave/northridge/nor_nwh090_a.txt','r+');
	fidv(4)=fopen('d:/berkleywave/northridge/nor_nwh090_v.txt','r+');
	ddt(4)=0.02;
	%next one 176
	fida(5)=fopen('d:/berkleywave/northridge/nor_nwh360_a.txt','r+');
	fidv(5)=fopen('d:/berkleywave/northridge/nor_nwh360_v.txt','r+');
	ddt(5)=0.02;
	%next one 215
	fida(6)=fopen('d:/berkleywave/northridge/nor_rrs228_a.txt','r+');
	fidv(6)=fopen('d:/berkleywave/northridge/nor_rrs228_v.txt','r+');
	ddt(6)=0.005;
	%next one 223
	fida(7)=fopen('d:/berkleywave/northridge/nor_sce018_a.txt','r+');
	fidv(7)=fopen('d:/berkleywave/northridge/nor_sce018_v.txt','r+');
	ddt(7)=0.005;
	%next one 255
	fida(8)=fopen('d:/berkleywave/northridge/nor_syl090_a.txt','r+');
	fidv(8)=fopen('d:/berkleywave/northridge/nor_syl090_v.txt','r+');
	ddt(8)=0.02;
	%next one 256
	fida(9)=fopen('d:/berkleywave/northridge/nor_syl360_a.txt','r+');
	fidv(9)=fopen('d:/berkleywave/northridge/nor_syl360_v.txt','r+');
	ddt(9)=0.02;
	%next one 289
	fida(10)=fopen('d:/berkleywave/northridge/nor_wpi046_a.txt','r+');
	fidv(10)=fopen('d:/berkleywave/northridge/nor_wpi046_v.txt','r+');
	ddt(10)=0.01;
	%=加載地震波=---------------------------------------------------------------
	nw=10;
    eb=0.333;
    um=zeros(1,sp);
    rm=zeros(1,sp);
	for i=1:nw
		dt=ddt(i); 
		ag=fscanf(fida(i),'%f',[1 inf])*9.81; 
		vg=fscanf(fidv(i),'%f',[1 inf])/100; 
		status=fclose(fida(i));status=fclose(fidv(i));
		st=length(ag);%
		at=dt*(st-1); %總時間為at；
		t=0:dt:at;
		cmag=max(abs(ag));
		%ag=ag*2.2/cmag;%%%為和肖結果比較，特添加；
		ag=ag*0.4*9.81/cmag;;%%%和CHAI比較，假定為0.4g，調整mvg；
		vg=vg*0.4*9.81/cmag;
		mag=max(abs(ag));mvg=max(abs(vg));
		acv=mag/9.81/mvg;
		tc=2*pi*2.0/2.5*mvg/mag;%計算場地特征周期Tc；cv=2.0;ca=2.5;
		%%%&***&***能量分析應采用相同震級地面運動加速度進行標準化；
		%%%計算強震時間td;td=t0.95-t0.05;
		ia=zeros(1,st);
		for m=2:st
            ia(m)=ia(m-1)+pi/2/9.81*(ag(m-1)^2+ag(m)^2)*dt/2;
		end
		iac1=0.1;
		iac2=0.1;
		for n=1:st
            ia1=abs(ia(n)/ia(st)-0.05);
            ia2=abs(ia(n)/ia(st)-0.95);
            if(ia1<iac1)
                iac1=ia1;
                t1=(n-1)*dt;
            end
            if(ia2<iac2)
                iac2=ia2;
                t2=(n-1)*dt;
            end
		end
		td=t2-t1;
		%%%做線彈性計算；
		para1=1;
		[rd1,ad1,rv1,av1,ra1,aa1]=response_l(para1,per,ag,dt,at);
		% figure(1)
		% mesh(per,t,ad1) %繪制出絕對位移時頻反應譜圖；
		% grid on
		% xlabel('周期 / 秒');
		% ylabel('時間 / 秒');
		[row,rank]=size(aa1);
		sa=zeros(1,rank);
		for ri=1:rank
            for rj=1:row
                if (abs(aa1(rj,ri))>sa(ri))
                    sa(ri)=aa1(rj,ri);%彈性加速度反應譜sa（最大絕對加速度反應）；
                end
            end
		end
		%%%做彈塑性計算；
		para1=1;
		[u1,nrd1,nad1,nrv1,nav1,nra1,naa1,rep]=response_n(para1,sa,per,ag,dt,at,eb);
		% figure(2)
		% mesh(per,t,nad1) %繪制出彈塑性位移時頻反應譜；
		% grid on
		% xlabel('周期 / 秒');
		% ylabel('時間 / 秒');
		um1=zeros(1,sp);
		for np=1:sp
            for nt=1:st
                if (abs(u1(nt,np))>um1(np))
                    um1(np)=abs(u1(nt,np));
                end
            end
		end
		
		um=um+um1;
		rm=rm+rep;
	end
	um=um/nw;
	rm=rm/nw;
    if(iy==1)
		figure(1)
        for icc=5:sp-1%%%誤差調整;
            if(abs(um(icc)-um(icc-1))>abs(um(icc-1)-um(icc-2))*3)
                um(icc)=(um(icc-1)+um(icc+1))/2;
            end
        end
		plot(per(2:sp),um(2:sp),'kd-.') %繪制出彈塑性位移延性反應譜；
		grid on
		xlabel('周期T(s)');
		ylabel('位移延性系數u');
		title('彈塑性位移延性反應譜');
        hold on;
		figure(2)
        for icc=5:sp-1%%%誤差調整;
            if(abs(rm(icc)-rm(icc-1))>abs(rm(icc-1)-rm(icc-2))*3)
                rm(icc)=(rm(icc-1)+rm(icc+1))/2;
            end
        end
		plot(per,rm,'kd-.') %繪制出累積滯回耗能譜；
		grid on
		xlabel('周期T(s)');
		ylabel('滯回耗能Eh(J/kg)');
		title('彈塑性滯回耗能譜');
        hold on;
    elseif(iy==2)
        figure(1)
        for icc=5:sp-1%%%誤差調整;
            if(abs(um(icc)-um(icc-1))>abs(um(icc-1)-um(icc-2))*3)
                um(icc)=(um(icc-1)+um(icc+1))/2;
            end
        end
		plot(per(2:sp),um(2:sp),'kh-.') %繪制出彈塑性位移延性反應譜；
        hold on;
        figure(2)
        for icc=5:sp-1%%%誤差調整;
            if(abs(rm(icc)-rm(icc-1))>abs(rm(icc-1)-rm(icc-2))*3)
                rm(icc)=(rm(icc-1)+rm(icc+1))/2;
            end
        end
		plot(per,rm,'kh-.') %繪制出累積滯回耗能譜；
        hold on;
    elseif(iy==3)
        figure(1)
		plot(per(2:sp),um(2:sp),'ko-.') %繪制出彈塑性位移延性反應譜；
        hold on;
        figure(2)
        for icc=5:sp-1%%%誤差調整;
            if(abs(rm(icc)-rm(icc-1))>abs(rm(icc-1)-rm(icc-2))*3)
                rm(icc)=(rm(icc-1)+rm(icc+1))/2;
            end
        end
		plot(per,rm,'ko-.') %繪制出累積滯回耗能譜；
        hold on;
    elseif(iy==4)
        figure(1)
        for icc=5:sp-1%%%誤差調整;
            if(abs(um(icc)-um(icc-1))>abs(um(icc-1)-um(icc-2))*3)
                um(icc)=(um(icc-1)+um(icc+1))/2;
            end
        end
		plot(per(2:sp),um(2:sp),'kp-.') %繪制出彈塑性位移延性反應譜；
        hold on;
        figure(2)
        for icc=5:sp-1%%%誤差調整;
            if(abs(rm(icc)-rm(icc-1))>abs(rm(icc-1)-rm(icc-2))*3)
                rm(icc)=(rm(icc-1)+rm(icc+1))/2;
            end
        end
		plot(per,rm,'kp-.') %繪制出累積滯回耗能譜；
        hold on;
    else
        figure(1)
        for icc=5:sp-1%%%誤差調整;
            if(abs(um(icc)-um(icc-1))>abs(um(icc-1)-um(icc-2))*3)
                um(icc)=(um(icc-1)+um(icc+1))/2;
            end
        end
		plot(per(2:sp),um(2:sp),'ks-.') %繪制出彈塑性位移延性反應譜；
        hold on;
        figure(2)
        for icc=5:sp-1%%%誤差調整;
            if(abs(rm(icc)-rm(icc-1))>abs(rm(icc-1)-rm(icc-2))*3)
                rm(icc)=(rm(icc-1)+rm(icc+1))/2;
            end
        end
		plot(per,rm,'ks-.') %繪制出累積滯回耗能譜；
        hold on;
    end       
end