inecal.m

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

function inecal

global inetextaxes;
global numero;
global A;
global Aine;
global inegraf;
global ninegraf;
global nninegraf;
global inedeltat;
global ineini;
global inedt;
global ineXi;
global inetf;
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;

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;

dt=deltat(inegraf);
k=0;
Nfft=2^(fix(log(Nn)/log(2))+1);

NNine=fix((inetf-10*dt)/inedt+1);
ineini=10*dt;
inetf=10*dt+(NNine-1)*inedt;
Q=[1.5 2 3 4 8 12 16];

for i=ineini:inedt:inetf;   
    k=k+1;
    f=1/i;
    w=2*pi/i;
    Bfft=fft(A(No+Nizq:Nf),Nfft);
    Tf1=1/w/w./((1.-([0:Nfft/2]/dt/Nfft/f).*([0:Nfft/2]/dt/Nfft/f))+((j*2*ineXi).*([0:Nfft/2]/dt/Nfft/f)));
    Tf2=1/w/w./((1.-([-Nfft/2+1:-1]/dt/Nfft/f).*([-Nfft/2+1:-1]/dt/Nfft/f))+((j*2*ineXi).*([-Nfft/2+1:-1]/dt/Nfft/f)));
    Tf=[Tf1 Tf2]';
    Bfft=Bfft.*Tf;
    Afft=ifft(Bfft); 
    m=1;
    R(m,k)=w*w*max(abs(Afft));

%    k=k+1;
%    w=2*pi/i;
%    Xt0=0;
%    Vt0=0;
%    Sd=0;
%    for j=No+Nizq:Nf-1;
%        A0=A(j);
%        a =(A(j+1)-A(j))/dt;
%        D=-a/w/w;
%        C=(-A0+2*ineXi*a/w)/w/w;
%        G1=Xt0-C;
%        G2=(Vt0+ineXi*w*G1+a/w/w)/w;
%        Xt0=exp(-ineXi*w*dt)*(G1*cos(w*dt)+G2*sin(w*dt))+C+D*dt;
%        Vt0=exp(-ineXi*w*dt)*(-G1*w*sin(w*dt)+G2*w*cos(w*dt))-ineXi*w*exp(-ineXi*w*dt)*(G1*cos(w*dt)+G2*sin(w*dt))+D;
%        if Sd<abs(Xt0);Sd=abs(Xt0);end;
%    end;
%    m=1;
%    R(m,k)=w*w*Sd;

    HE(m,k)=0;
    Umax(m,k)=1;
    Xmax(m,k)=R(m,k)/w/w;
    for q=[Q];
        m=m+1;
        R(m,k)=R(1,k)/q;
        HE(m,k)=0;
        Umax(m,k)=0;
        X0=0;
        V0=0;
        A0=0;
        R0=0;
        P0=0;
        tramo=1;
        for l=No+Nizq:Nf-1;
            if tramo==1;
               DP=(A(l)-P0)-A0*dt/2*(4*ineXi*w+dt*w*w)-w*w*dt*V0;
               DK=1+ineXi*w*dt+w*w*dt*dt/4;
               DA=DP/DK;
               Xu=X0+dt*dt/4*(DA+2*A0)+dt*V0;
               Ru=(Xu-X0)*w*w+R0;
               if abs(Ru)>=R(m,k);
                  caso=12;
                  Ru=Ru/abs(Ru)*R(m,k);
                  DX=(Ru-R0)/w/w;
                  c0=4*DX;
                  c1=4*(ineXi*w*DX-V0);
                  c2=w*w*DX-4*ineXi*w*V0-2*A0;
                  c3=(P0-A(l))/dt;
                  raices=roots([c3 c2 c1 c0]);
                  iraiz1=imag(raices(1));
                  rraiz1=real(raices(1)); 
                  if iraiz1==0 & rraiz1>=0 & rraiz1<=dt;
                     dtt=rraiz1;
                  else;
                     iraiz2=imag(raices(2));
                     rraiz2=real(raices(2));
                     if iraiz2==0 & rraiz2>=0 & rraiz2<=dt;
                        dtt=rraiz2;
                     else;
                        rraiz3=real(raices(3));
                        dtt=rraiz3;
                     end;
                  end;
                  DP=(A(l)-P0)/dt*dtt-A0*dtt/2*(4*ineXi*w+dtt*w*w)-w*w*dtt*V0;
                  DK=1+ineXi*w*dtt+w*w*dtt*dtt/4;
                  DA=DP/DK;
                  Xu=X0+dtt*dtt/4*(DA+2*A0)+dtt*V0;
                  HE(m,k)=HE(m,k)+(Ru+R0)*(Xu-X0)/2;
                  R0=Ru;
                  X0=Xu;
                  V0=V0+dtt/2*(DA+2*A0);
                  A0=A0+DA;             
                  P0=P0+(A(l)-P0)/dt*dtt;
                  dtt=dt-dtt;
                  DA=A(l)-P0;
                  Xu=X0+dtt*dtt/4*(DA+2*A0)+dtt*V0;
                  Vu=V0+dtt/2*(DA+2*A0);
                  if Vu*V0<=0;
                     caso=21;
                     DV=-V0;
                     c0=2*DV;
                     c1=-2*A0;
                     c2=(P0-A(l))/dtt;
                     raices=roots([c2 c1 c0]);
                     iraiz1=imag(raices(1));
                     rraiz1=real(raices(1)); 
                     if iraiz1==0 & rraiz1>=0 & rraiz1<=dtt;
                        dttt=rraiz1;
                     else
                        rraiz2=real(raices(2));
                        dttt=rraiz2;
                     end;
                     DA=(A(l)-P0)/dtt*dttt;
                     Xu=X0+dttt*dttt/4*(DA+2*A0)+dttt*V0;
                     if abs(Xu)/R(m,k)*w*w>Umax(m,k);
                        Umax(m,k)=abs(Xu)/R(m,k)*w*w;
                        Xmax(m,k)=abs(Xu);
                     end;
                     HE(m,k)=HE(m,k)+R0*(Xu-X0);
                     X0=Xu;
                     V0=0;
                     A0=A0+DA;             
                     P0=P0+(A(l)-P0)/dtt*dttt;
                     dttt=dtt-dttt;
                     DP=(A(l)-P0)-A0*dttt/2*(4*ineXi*w+dttt*w*w)-w*w*dttt*V0;
                     DK=1+ineXi*w*dttt+w*w*dttt*dttt/4;
                     DA=DP/DK;
                     Xu=X0+dttt*dttt/4*(DA+2*A0)+dttt*V0;
                     Ru=(Xu-X0)*w*w+R0;
                     if abs(Ru)>=R(m,k);
                        axes(inetextaxes);hold off;
                        axis('off')
                        text('Position',[0.05 0.05],'Units','normalized','Color','r','String',['ERROR T=' num2str(i) ' Q=' num2str(q)]);
                        axes(ineaxes);
                        hold off;
                     end;   
                     HE(m,k)=HE(m,k)+(Ru+R0)*(Xu-X0)/2;
                     R0=Ru;
                     X0=Xu;
                     V0=V0+dttt/2*(DA+2*A0);
                     A0=A0+DA;             
                     P0=A(l);
                     tramo=1;                              
                  else;
                     HE(m,k)=HE(m,k)+R0*(Xu-X0);
                     X0=Xu;
                     V0=V0+dtt/2*(DA+2*A0);
                     A0=A0+DA;             
                     P0=A(l);
                     tramo=2;                           
                  end;
               else;
                  HE(m,k)=HE(m,k)+(Ru+R0)*(Xu-X0)/2;
                  R0=Ru;
                  X0=Xu;
                  V0=V0+dt/2*(DA+2*A0);
                  A0=A0+DA;             
                  P0=A(l);                  
               end;                   
            else;
               DA=A(l)-P0;
               Xu=X0+dt*dt/4*(DA+2*A0)+dt*V0;
               Vu=V0+dt/2*(DA+2*A0);
               if Vu*V0<=0;
                  caso=21;
                  DV=-V0;
                  c0=2*DV;
                  c1=-2*A0;
                  c2=(P0-A(l))/dt;
                  raices=roots([c2 c1 c0]);
                  iraiz1=imag(raices(1));
                  rraiz1=real(raices(1)); 
                  if iraiz1==0 & rraiz1>=0 & rraiz1<=dt;
                     dtt=rraiz1;
                  else
                     rraiz2=real(raices(2));
                     dtt=rraiz2;
                  end;
                  DA=(A(l)-P0)/dt*dtt;
                  Xu=X0+dtt*dtt/4*(DA+2*A0)+dtt*V0;
                  if abs(Xu)/R(m,k)*w*w>Umax(m,k);
                     Umax(m,k)=abs(Xu)/R(m,k)*w*w;
                     Xmax(m,k)=abs(Xu);
                  end;
                  HE(m,k)=HE(m,k)+R0*(Xu-X0);
                  X0=Xu;
                  V0=0;
                  A0=A0+DA;             
                  P0=P0+(A(l)-P0)/dt*dtt;
                  dtt=dt-dtt;
                  DP=(A(l)-P0)-A0*dtt/2*(4*ineXi*w+dtt*w*w)-w*w*dtt*V0;
                  DK=1+ineXi*w*dtt+w*w*dtt*dtt/4;
                  DA=DP/DK;
                  Xu=X0+dtt*dtt/4*(DA+2*A0)+dtt*V0;
                  Ru=(Xu-X0)*w*w+R0;
                  if abs(Ru)>=R(m,k);
                     axes(inetextaxes);hold off;
                     axis('off')
                     text('Position',[0.05 0.05],'Units','normalized','Color','r','String',['ERROR T=' num2str(i) ' Q=' num2str(q)]);
                     axes(ineaxes);
                     hold off;
                  end;   
                  HE(m,k)=HE(m,k)+(Ru+R0)*(Xu-X0)/2;
                  R0=Ru;
                  X0=Xu;
                  V0=V0+dtt/2*(DA+2*A0);
                  A0=A0+DA;             
                  P0=A(l);
                  tramo=1;                              
               else;
                  HE(m,k)=HE(m,k)+R0*(Xu-X0);
                  X0=Xu;
                  V0=V0+dt/2*(DA+2*A0);
                  A0=A0+DA;             
                  P0=A(l);                  
               end;                   
            end;      
        end;
        if Umax(m,k)>4 & q>=4;
           n=0;
           for l=m+1:8;
               n=n+1;
               Umax(l,k)=max([Umax(m,k)+0.1*n,Q(l-1)]);
               Xmax(l,k)=Xmax(m,k);
               HE(l,k)=HE(m,k);
               R(l,k)=R(1,k)/Q(l-1);
           end;
           break;
        end;
    end;
end;      

[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 1.5 2 3 4];

k=0;
for i=ineini:inedt:inetf;
    k=k+1;
    Sainel=[Sainel;(interp1(UUmax(:,k),RR(:,k),X))'];
    SHinel=[SHinel;(interp1(UUmax(:,k),HHE(:,k),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=[0 0.10 0.20 0.30 0.40];

k=0;
for i=ineini:inedt:inetf;
    k=k+1;
    IDPA  =[IDPA  ;(interp1(DDPA(:,k),RRR(:,k),XPA))'];
end;

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

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;

saingra;