inecal2.m

德国人开发的地震处理分析软件

function inecal2

global inetextaxes;
global numero;
global A;
global Aine;
global inegraf;
global ninegraf;
global nninegraf;
global inedeltat;
global ineXi;
global ineBeta;
global ineUmon;
global deltat;
global Limiteizq;
global Limiteder;
global numeroine;
global Umax;
global UUmax;
global R;
global RR;
global HE;
global HHE;
global Xmax;
global XXmax;
global NNine;
global inedt;
global directorio;
global directorio1;

nninegraf=inegraf;

No=numero(inegraf);
Nn=numero(inegraf+1)-No;

if Limiteizq==Limiteder;
   Nizq=0;
   Nder=Nn;
else;
   Nizq=fix(Limiteizq/deltat(inegraf));
   Nder=fix(Limiteder/deltat(inegraf));
end;

if Nizq>Nder;Nizq=Mud;Nizq=Nder;Nder=Mud;end;

if Nder>=Nn;
   Nf=No+Nn-1;
else;
   Nf=No+Nder;
   Nn=Nder-Nizq+1;
end;

%cd d:\matlab\toolbox\prima
cd (directorio1);

fid=fopen('Acel.dat','wt');
fprintf(fid,'%12.4e\n',A(No+Nizq:Nf));
fclose(fid);

Nacel=Nf-No-Nizq+1;
dt=deltat(inegraf);
ineini=10*dt;
inetf=10*dt+(NNine-1)*inedt;

% Si se cambia el n鷐ero de resistencias a calcular, NQ, debe cambiarse 
% no s髄o el vector Q, sino tambi閚 el formato de impresi髇 del mismo.

% El mismo comentario es v醠ido si se cambia el n鷐ero de ductilidades
% a calcular, ND.

NQ=12;
Q=[2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 8.0 12.0 18.0];

ND=4;
DUC=[1.5 2.0 3.0 4.0];

% Si se cambia el n鷐ero de ductilidades a calcular debe cambiarse
% tambi閚 el n鷐ero de 韓dices de Park & Ang a calcular y el formato
% para el c醠culo de los vector RUmax y RIDPA.

XPA=[0.10 0.20 0.30 0.40];
XPA=[0.0 XPA];

fid=fopen('Sdat.dat','wt');
fprintf(fid,'%8i%8.3f%8i%8i%8i%8.3f%8.3f%8.3f\n',Nacel,dt,NQ,ND,NNine,inedt,10*dt,ineXi);
fprintf(fid,'%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f%8.3f\n',Q'); 
fprintf(fid,'%8.3f%8.3f%8.3f%8.3f\n',DUC'); 
fclose(fid);

!inecal

load Umax.sal;
load HE.sal;
load R.sal;
load Xmax.sal;

Umax=Umax';
HE=HE';
R=R';
Xmax=Xmax';

delete Acel.dat
delete Sdat.dat
delete Umax.sal
delete HE.sal
delete R.sal
delete Xmax.sal

[UUmax,I]=sort(Umax);

k=0;
for i=ineini:inedt:inetf;
    k=k+1;
    l=0;
    for j=[I(:,k)'];
        l=l+1;
        RR(l,k)=R(j,k);
        HHE(l,k)=HE(j,k);
        XXmax(l,k)=Xmax(j,k);
    end;
end;           

X=[1.0 DUC];

Dmin=min(UUmax(length(Q)+1,:));

for i=1:length(X);
    if Dmin<X(i);
       XX=X(i-1)*ones(1,length(X)-i+1);
       X=[X(1:i-1) XX];
       i=length(X)+1;
    end;
end;
       
for i=1:NNine;
    Sainel=[Sainel;(interp1(UUmax(:,i),RR(:,i),X))'];
    SHinel=[SHinel;(interp1(UUmax(:,i),HHE(:,i),X))'];
end;

k=0;
for i=ineini:inedt:inetf;
    k=k+1;
    DPA  =[DPA  ,(UUmax(:,k)/ineUmon+ineBeta*HHE(:,k)./RR(:,k)/ineUmon./XXmax(:,k).*UUmax(:,k))];
end;

DPA(1,:)=DPA(1,:)*0;

[DDPA,I]=sort(DPA);

k=0;
for i=ineini:inedt:inetf;
    k=k+1;
    l=0;
    for j=[I(:,k)'];
        l=l+1;
        RRR(l,k)=RR(j,k);
    end;
end;

XPA(1)=DDPA(1,1);

for i=1:NNine;
    IDPA  =[IDPA  ;(interp1(DDPA(:,i),RRR(:,i),XPA))'];
end;

Sdinel=Xmax(1:5,:)';


% cambiar cuando se cambie el n鷐ero de ductilidades a calcular


RUmax=[R(1,:)' R(1,:)' R(1,:)' R(1,:)' R(1,:)']./Sainel;
RIDPA=[R(1,:)' R(1,:)' R(1,:)' R(1,:)' R(1,:)']./IDPA;

ind=0;
j=0;
for i=ninegraf;
    j=j+1;
    if i==inegraf;
       inedeltat(j)=inedt;
       Aine=[Aine(1:numeroine(j)-1,:);[Sainel;Sdinel;SHinel;RUmax;IDPA;RIDPA];...
             Aine(numeroine(j+1):length(Aine),:)];
       dif=6*NNine-numeroine(j+1)+numeroine(j);
       if j==length(ninegraf);
          numeroine=[numeroine(1:j);numeroine(j)+6*NNine];
       else; 
          numeroine=[numeroine(1:j);numeroine(j)+6*NNine;...
                    (numeroine(j+2:length(numeroine))+dif)];
       end;
       ind=1;
       break; 
    end;
end;

if ind==0;
   inedeltat=[inedeltat;inedt];
   Aine=[Aine;Sainel;Sdinel;SHinel;RUmax;IDPA;RIDPA];
   ninegraf=[ninegraf inegraf];
   numeroine=[numeroine;numeroine(length(numeroine))+6*NNine];
end;

cd (directorio);

saingra;