rescal.m

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

function rescal

global numero;
global A;
global Ares;
global resgraf;
global nresgraf;
global nnresgraf;
global resdeltat;
global resdt;
global resXi;
global restf;
global deltat;
global Limiteizq;
global Limiteder;
global numerores;
global directorio;
global directorio1;

nnresgraf=resgraf;

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

if Limiteizq==Limiteder;
   Nizq=0;
   Nder=Nn;
else
   Nizq=fix(Limiteizq/deltat(resgraf));
   Nder=fix(Limiteder/deltat(resgraf));
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 (directorio1)

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

Nacel=Nf-No-Nizq+1;
dt=deltat(resgraf);
NNres=fix(restf/resdt);

fid=fopen('Sdat.dat','wt');
fprintf(fid,'%8i%8.3f%8i%8.3f%8.3f%8.3f',Nacel,dt,NNres,resdt,resdt,resXi);
fclose(fid);

!ecal

load ASd.sal;
load ASv.sal;
load ASa.sal;

delete Acel.dat
delete Sdat.dat
delete ASd.sal
delete ASv.sal
delete ASa.sal

%k=0;
%for i=resdt:resdt:restf;   
%    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*resXi*a/w)/w/w;
%        G1=Xt0-C;
%        G2=(Vt0+resXi*w*G1+a/w/w)/w;
%        Xt0=exp(-resXi*w*DT)*(G1*cos(w*DT)+G2*sin(w*DT))+C+D*DT;
%        Vt0=exp(-resXi*w*DT)*(-G1*w*sin(w*DT)+G2*w*cos(w*DT))-resXi*w*exp(-resXi*w*DT)*(G1*cos(w*DT)+G2*sin(w*DT))+D;
%        if Sd<abs(Xt0);Sd=abs(Xt0);end;
%    end;
%    ASd(k)=Sd;        
%    ASv(k)=Sd*w;
%    ASa(k)=Sd*w*w;
%end;


%k=0;
%Nfft=2^(fix(log(Nn)/log(2))+1);
%for l=resdt:resdt:restf;   
%    k=k+1;
%    f=1/l;
%    w=2*pi*f;
%    Bfft=fft(A(No+Nizq:Nf),Nfft);
%    Tf1=1/w/w./((1.-([0:Nfft/2]/DT/Nfft/f).*([0:Nfft/2]/DT/Nfft/f))+((i*2*resXi).*([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))+((i*2*resXi).*([-Nfft/2+1:-1]/DT/Nfft/f)));
%    Tf=[Tf1 Tf2]';
%    Bfft=Bfft.*Tf;
%    Afft=ifft(Bfft);
%    Xmax=max(abs(Afft));
%    ASd(k)=Xmax;      
%    ASv(k)=ASd(k)*w;
%    ASa(k)=ASd(k)*w*w;
%end;

ind=0;
j=0;
for i=nresgraf;
    j=j+1;
    if i==resgraf;
       resdeltat(j)=resdt;
       Ares=[Ares(1:numerores(j)-1,:);[ASd ASv ASa];...
             Ares(numerores(j+1):length(Ares),:)];
       dif=NNres-numerores(j+1)+numerores(j);
       if j==length(nresgraf);
          numerores=[numerores(1:j);numerores(j)+NNres];
       else; 
          numerores=[numerores(1:j);numerores(j)+NNres;...
                    (numerores(j+2:length(numerores))+dif)];
       end;
       ind=1;
       break; 
    end;
end;

if ind==0;
   resdeltat=[resdeltat;resdt];
   Ares=[Ares;[ASd ASv ASa]];
   nresgraf=[nresgraf resgraf];
   numerores=[numerores;numerores(length(numerores))+NNres];
end;

cd (directorio);

sagra;