% SLB Project 
% Program to investigate the time dependence 
% of rotation rate energy 

% Pham D. N., H. Igel, and M. Kaeser, 
% Geophysics, LMU Munich.
% Jan. 2009


clear all;
close all;

mtype=input('Enter the type of the medium: 1 or 2 (1- Homo, 2 - Layer)):');
        if mtype==1
            f_name1='Homo';
            Res=input('Enter the source index (1 to 5):');
        else
            f_name1='Layer';
            Res=input('Enter the source index (1 or 5):');
        end

        f_name2=[num2str(Res)];

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Scale factor for the source so that in all case Mw=0; 
        if (Res==1 & mtype==1)
            Amp_coef=8.4074865549767914e+005; 
        else
            Amp_coef=1;
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
N_Rec=input('Enter the receiver index (1 to 15):');
f_name3=[num2str(N_Rec)];

        if Res<10
            filename=[f_name1 '0' f_name2 ' (' f_name3 ').dat'];
        else
            filename=[f_name1 f_name2 ' (' f_name3 ').dat'];
        end


        




cd ..; % Shift to data directory


numvars=6; % Number of components;

% Load the data of the simulated seismograms
fid = fopen(filename,'r');
tmp=fgets(fid);
tmp=fgets(fid);
tmp=fgets(fid);
tmp=fgets(fid);
tmp=fgets(fid);
S=fscanf(fid,'%f');
fclose('all');

for j=1:length(S)/(numvars+1)
    for i=1:numvars+1
         SEIS(i,j)=S(i+(j-1)*(numvars+1));
    end
end
SEIS(2:i,:)=SEIS(2:i,:)*Amp_coef;


    
Ttime=SEIS(1,:);
Vn=SEIS(2,:);
Ve=SEIS(3,:);
Vz=SEIS(4,:);
Rn=SEIS(5,:);
Re=SEIS(6,:);
Rz=SEIS(7,:);

dts=Ttime(2)-Ttime(1);
nts=length(Ttime);


    Ae= diff(Ve)/dts;Ae(length(Ae)+1)=0;
    An= diff(Vn)/dts;An(length(An)+1)=0;
    Az= diff(Vz)/dts;Az(length(Az)+1)=0;
    

sig(1,:)=Vn;
sig(2,:)=Ve;
sig(3,:)=Vz;

sig(4,:)=An;
sig(5,:)=Ae;
sig(6,:)=Az;

sig(7,:)=Rn;
sig(8,:)=Re;
sig(9,:)=Rz;



noise_amp_P=0; % N/S in percentage;
noise_amp=noise_amp_P/100; % N/S in decimal;

for nnn=1:numvars+3
    % Generate a series of random noise 
    pert(nnn,:)=(rand(1,length(Ttime))-0.5).*2;
    pert(nnn,:)=noise_amp*max(abs(sig(nnn,:))).*pert(nnn,:);
    % Add to the signals
    sign(nnn,:)=sig(nnn,:)+pert(nnn,:);
end


sign(numvars+4,:)=  (sign(4,:).^2+sign(5,:).^2+sign(6,:).^2); % E_A   
sign(numvars+5,:)=  (sign(7,:).^2+sign(8,:).^2+sign(9,:).^2); % E_R   

                
        VARS=['V_N';'V_E';'V_Z';...
        'A_N';'A_E';'A_Z';...
        'R_N';'R_E';'R_Z';...
        'Ea ';'Er '];

    
    
k=[1:3 7:11];
    
    for ii=1:length(k);
        i=k(ii);

            if (ii ==1 | ii==4)
                ms='k';
            elseif (ii ==2 | ii==5 | ii==8)
                ms='b';
            else
                ms='r';
            end
        
                ghy(i)=max(abs(sign(i,:)));
        
        if ii<= 3   
            subplot(4, 1, 1), 
            hold on; plot(Ttime,sign(i,:)/ghy(i)+(ii-1),ms);set(gca,'YTick',[]);
            text(0.02, (ii-1)+0.4,VARS(i,:));
            axis([min(Ttime) max(Ttime) -1 3]); box on;
        elseif (ii>3 & ii<=6)
            subplot(4, 1, 2), 
            hold on, plot(Ttime,sign(i,:)/ghy(i)+(ii-4),ms);set(gca,'YTick',[]);
            text(0.02, (ii-4)+0.4,VARS(i,:));
            axis([min(Ttime) max(Ttime) -1 3]); box on;
        else
            subplot(4, 1, ii-4),
            plot(Ttime, sign(i,:)/ghy(i),ms),set(gca,'YTick',[]);
            axis([min(Ttime) max(Ttime) -0.1 1]); box on;
            text(0.02, 0.1,VARS(i,:));

        end

        
       
    end
    xlabel('Time (s)');
    
%pfile=sprintf('./Figs/Fig6_SPSeparate_%i_%i_%i_%i.jpg',mtype,Res,N_Rec,noise_amp_P);eval([' print -djpeg100 ' pfile ]);


cd MatLabCodes; % Back to program directory

