ПРИЛОЖЕНИЕ Б

(обязательное)

Построение волновых спектрограмм землетрясений

Программа по построению волновых форм спектрограмм землетрясений и микросейсм до и после землетрясений

p='f:\ATR\';

p2='f:\';

%nomer='13';

switchnomer

case '01'

filenameX1='2014204214130320_POLE__1_2.atr';

filenameY1='2014204214130320_POLE__1_3.atr';

filenameX2='2014204224130320_POLE__1_2.atr';

filenameY2='2014204224130320_POLE__1_3.atr';

filetime=datenum(2014,07,23,21,41,30.320);

begtime=datenum(2014,07,23,22,48,17);

ZTname='M=1.1,Az=000';

timeDo=50; timePosle=100;

z1=-125; z2=-90;

case '02'

filenameX1='2014205135635320_POLE__1_2.atr';

filenameY1='2014205135635320_POLE__1_3.atr';

filenameX2='2014205145635320_POLE__1_2.atr';

filenameY2='2014205145635320_POLE__1_3.atr';

filetime=datenum(2014,07,24,13,56,35.320);

begtime=datenum(2014,07,24,14,48,9); %!!!!!

ZTname='M=2.4,Az=097';

timeDo=50; timePosle=100;

z1=-120; z2=-80;

Продолжение приложения Б

filenameX1='2014205165635320_POLE__1_2.atr';

filenameY1='2014205165635320_POLE__1_3.atr';

filenameX2='2014205175635320_POLE__1_2.atr';

filenameY2='2014205175635320_POLE__1_3.atr';

filetime=datenum(2014,07,24,16,56,35.320);

begtime=datenum(2014,07,24,17,12,3);

ZTname='M=3.8,Az=027';

timeDo=50; timePosle=250;

z1=-120; z2=-50;

case '05'

filenameX1='2014205175635320_POLE__1_2.atr';

filenameY1='2014205175635320_POLE__1_3.atr';

filenameX2='2014205185635320_POLE__1_2.atr';

filenameY2='2014205185635320_POLE__1_3.atr';

filetime=datenum(2014,07,24,17,56,35.320);

begtime=datenum(2014,07,24,18,40,42);

ZTname='M=1.8,Az=107';

timeDo=50; timePosle=100;

z1=-125; z2=-80;

case '06'

filenameX1='2014206055635320_POLE__1_2.atr';

filenameY1='2014206055635320_POLE__1_3.atr';

filenameX2='2014206065635320_POLE__1_2.atr';

filenameY2='2014206065635320_POLE__1_3.atr';

filetime=datenum(2014,07,25,05,56,35.320);

begtime=datenum(2014,07,25,06,07,2);

ZTname='M=0.9,Az=239';

timeDo=50; timePosle=70;

z1=-115; z2=-95;

Продолжение приложения Б

case '07'

filenameX1='2014206115635320_POLE__1_2.atr';

filenameY1='2014206115635320_POLE__1_3.atr';

filenameX2='2014206125635320_POLE__1_2.atr';

filenameY2='2014206125635320_POLE__1_3.atr';

filetime=datenum(2014,07,25,11,56,35.320);

begtime=datenum(2014,07,25,12,43,6);

ZTname='M=0.7,Az=229';

timeDo=50; timePosle=40;

z1=-125; z2=-80;

case '08'

filenameX1='2014206175635320_POLE__1_2.atr';

filenameY1='2014206175635320_POLE__1_3.atr';

filenameX2='2014206185635320_POLE__1_2.atr';

filenameY2='2014206185635320_POLE__1_3.atr';

filetime=datenum(2014,07,25,17,56,35.320);

begtime=datenum(2014,07,25,18,04,27);

ZTname='M=2.0,Az=117';

timeDo=50; timePosle=150;

z1=-125; z2=-85;

case '09'

filenameX1='2014206205635320_POLE__1_2.atr';

filenameY1='2014206205635320_POLE__1_3.atr';

filenameX2='2014206215635320_POLE__1_2.atr';

filenameY2='2014206215635320_POLE__1_3.atr';

filetime=datenum(2014,07,25,20,56,35.320);

begtime=datenum(2014,07,25,21,24,7);

ZTname='M=1.6,Az=249';

timeDo=50; timePosle=100;

Продолжение приложения Б

z1=-125; z2=-95;

case '10'

filenameX1='2014207175635320_POLE__1_2.atr';

filenameY1='2014207175635320_POLE__1_3.atr';

filenameX2='2014207185635320_POLE__1_2.atr';

filenameY2='2014207185635320_POLE__1_3.atr';

filetime=datenum(2014,07,26,17,56,35.320);

begtime=datenum(2014,07,26,18,06,31);

ZTname='M=1.7,Az=027';

timeDo=50; timePosle=100;

z1=-125; z2=-85;

case '11'

filenameX1='2014208055635320_POLE__1_2.atr';

filenameY1='2014208055635320_POLE__1_3.atr';

filenameX2='2014208065635320_POLE__1_2.atr';

filenameY2='2014208065635320_POLE__1_3.atr';

filetime=datenum(2014,07,27,05,56,35.320);

begtime=datenum(2014,07,27,06,32,58);

ZTname='M=2.1,Az=237';

timeDo=50; timePosle=130;

z1=-120; z2=-90;

case '12'

filenameX1='2014208075635320_POLE__1_2.atr';

filenameY1='2014208075635320_POLE__1_3.atr';

filenameX2='2014208085635320_POLE__1_2.atr';

filenameY2='2014208085635320_POLE__1_3.atr';

filetime=datenum(2014,07,27,07,56,35.320);

begtime=datenum(2014,07,27,08,24,25);

ZTname='M=1.7,Az=217';

Продолжение приложения Б

timeDo=50; timePosle=115;

z1=-125; z2=-95;

case '13'

filenameX1='2014208125635320_POLE__1_2.atr';

filenameY1='2014208125635320_POLE__1_3.atr';

filenameX2='2014208135635320_POLE__1_2.atr';

filenameY2='2014208135635320_POLE__1_3.atr';

filetime=datenum(2014,07,27,12,56,35.320);

begtime=datenum(2014,07,27,13,12,9);

ZTname='M=4.0,Az=329';

timeDo=50; timePosle=300;

z1=-125; z2=-55;

case '14'

filenameX1='2014210225635320_POLE__1_2.atr';

filenameY1='2014210225635320_POLE__1_3.atr';

filenameX2='2014210235635320_POLE__1_2.atr';

filenameY2='2014210235635320_POLE__1_3.atr';

filetime=datenum(2014,07,29,22,56,35.320);

begtime=datenum(2014,07,29,23,32,35);

ZTname='M=2.1,Az=269';

timeDo=50; timePosle=150;

z1=-125; z2=-90;

case '15'

filenameX1='2014210225635320_POLE__1_2.atr';

filenameY1='2014210225635320_POLE__1_3.atr';

filenameX2='2014210235635320_POLE__1_2.atr';

filenameY2='2014210235635320_POLE__1_3.atr';

filetime=datenum(2014,07,29,22,56,35.320);

begtime=datenum(2014,07,29,23,35,26);

Продолжение приложения Б

ZTname='M=1.9,Az=247';

timeDo=50; timePosle=130;

z1=-125; z2=-95;

case '16'

filenameX1='2014211005635320_POLE__1_2.atr';

filenameY1='2014211005635320_POLE__1_3.atr';

filenameX2='2014211015635320_POLE__1_2.atr';

filenameY2='2014211015635320_POLE__1_3.atr';

filetime=datenum(2014,07,30,00,56,35.320);

begtime=datenum(2014,07,30,01,47,28);

ZTname='M=1.8,Az=217';

timeDo=50; timePosle=120;

z1=-125; z2=-85;

case '17'

filenameX1='2014211075635320_POLE__1_2.atr';

filenameY1='2014211075635320_POLE__1_3.atr';

filenameX2='2014211085635320_POLE__1_2.atr';

filenameY2='2014211085635320_POLE__1_3.atr';

filetime=datenum(2014,07,30,07,56,35.320);

begtime=datenum(2014,07,30,08,20,48);

ZTname='M=1.8,Az=099';

timeDo=50; timePosle=130;

z1=-125; z2=-80;

case '18'

filenameX1='2014212005635320_POLE__1_2.atr';

filenameY1='2014212005635320_POLE__1_3.atr';

filenameX2='2014212015635320_POLE__1_2.atr';

filenameY2='2014212015635320_POLE__1_3.atr';

filetime=datenum(2014,07,31,00,56,35.320);

Продолжение приложения Б

begtime=datenum(2014,07,31,01,54,13);

ZTname='M=1.7,Az=237';

timeDo=50; timePosle=110;

z1=-125; z2=-90;

case '19'

filenameX1='2014212025635320_POLE__1_2.atr';

filenameY1='2014212025635320_POLE__1_3.atr';

filenameX2='2014212035635320_POLE__1_2.atr';

filenameY2='2014212035635320_POLE__1_3.atr';

filetime=datenum(2014,07,31,02,56,35.320);

begtime=datenum(2014,07,31,03,43,5);

ZTname='M=2.7,Az=217';

timeDo=50; timePosle=120;

z1=-125; z2=-90;

case '20'

filenameX1='2014213035956320_POLE__1_2.atr';

filenameY1='2014213035956320_POLE__1_3.atr';

filenameX2='2014213045956320_POLE__1_2.atr';

filenameY2='2014213045956320_POLE__1_3.atr';

filetime=datenum(2014,08,01,03,59,56.320);

begtime=datenum(2014,08,01,04,46,0);

ZTname='M=2.1,Az=262';

timeDo=50; timePosle=150;

z1=-125; z2=-75;

case '21'

filenameX1='2014213035956320_POLE__1_2.atr';

filenameY1='2014213035956320_POLE__1_3.atr';

filenameX2='2014213045956320_POLE__1_2.atr';

filenameY2='2014213045956320_POLE__1_3.atr';

Продолжение приложения Б

filetime=datenum(2014,08,01,03,59,56.320);

begtime=datenum(2014,08,01,04,49,29);

ZTname='M=1.4,Az=266';

timeDo=30; timePosle=70;

z1=-125; z2=-90;

case '22'

filenameX1='2014213035956320_POLE__1_2.atr';

filenameY1='2014213035956320_POLE__1_3.atr';

filenameX2='2014213045956320_POLE__1_2.atr';

filenameY2='2014213045956320_POLE__1_3.atr';

filetime=datenum(2014,08,01,03,59,56.320);

begtime=datenum(2014,08,01,04,51,32);

ZTname='M=0.8,Az=262';

timeDo=30; timePosle=70;

z1=-125; z2=-100;

end

fmin=.3; dt=3;

delimeter=' ';

nStrok=9; L=86400; fs=200;

m=1.583e-6/2001.67*1000;

% load

data666 = importdata([p filenameX1], delimeter, nStrok);

dataX=data666.data;

data666 = importdata([p filenameX2], delimeter, nStrok);

dataX=[dataX; data666.data]*m;

data666 = importdata([p filenameY1], delimeter, nStrok);

dataY=data666.data;

data666 = importdata([p filenameY2], delimeter, nStrok);

dataY=[dataY; data666.data]*m;

Продолжение приложения Б

cleardata666;

% filter

[b,a]=butter(3,fmin/fs*2,'high');

dataXf=filter(b,a,dataX-mean(dataX));

dataYf=filter(b,a,dataY-mean(dataY));

% cut

t1=round((begtime-timeDo/L-filetime)*L*fs);

t2=round((begtime+timePosle/L-filetime)*L*fs);

X=dataXf(t1:t2); Y=dataYf(t1:t2);

% spectrogram

sp01(X,fs,1000,900,0,z1,z2)

shadinginterp, ylim([1,30]);

print(gcf,'-dpng',[p2 nomer '_' ZTname '_X.png'],'-r300')

close(gcf);

sp01(Y,fs,1000,900,0,z1,z2)

shadinginterp, ylim([1,30]);

print(gcf,'-dpng',[p2 nomer '_' ZTname '_Y.png'],'-r300')

close(gcf);

function sp01(fname,fs,wl,overlap, meaning,zmin,zmax)

%load file

ifischar(fname)

delimeter=' ';

nStrok=9;

data666 = importdata(fname, delimeter, nStrok);

data1=data666.data;

clear data666;

else

data1=fname;

end

Продолжение приложения Б

% computing

if meaning>1

data1=AntiTrendFast(data1,meaning);

end

[~,F,T,P]=spectrogram(data1,hann(wl),overlap,wl,fs);

F(1:3,:)=[];

P(1:3,:)=[];

%plot 1st

figure;

subplot('position',[0.04 0.75 0.94 0.22]);

set(gca,'fontSize',9)

plot((1:length(data1))./fs,data1); axis tight; grid on;

aaa=get(gca,'ylim');

aaa(1)=aaa(1)-0.02*(aaa(2)-aaa(1));

aaa(2)=aaa(2)+0.02*(aaa(2)-aaa(1));

Ylim(aaa);

%plot 2nd

subplot('position',[0.04 0.05 0.94 0.62]);

set(gca,'fontSize',9)

surf(T,F,10*log10(P),'edgecolor','none');

aa1=0; aa2=length(data1)/fs; aa3=F(1,1); aa22=F(size(F)); aa4=aa22(1);

axis ([aa1 aa2 aa3 aa4]);

set(gca,'yscale','log');

colorbar('east');

colormap(jet(4096));

set(gca,'clim',[zminzmax]);

holdon;

Продолжение приложения Б

a=[.0001 .0001 .0002 .0002 .0003 .0003 .0004 .0004 .0005 .0005 .0006 .0006 .0007 .0007 .0008 .0008 .0009 .0009 ...

.001 .001 .002 .002 .003 .003 .004 .004 .005 .005 .006 .006 .007 .007 .008 .008 .009 .009 ...

.01 .01 .02 .02 .03 .03 .04 .04 .05 .05 .06 .06 .07 .07 .08 .08 .09 .09 ...%9*2*6=108

.1 .1 .2 .2 .3 .3 .4 .4 .5 .5 .6 .6 .7 .7 .8 .8 .9 .9 ...

1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 9 ...

10 10 20 20 30 30 40 40 50 50 60 60 70 70 80 80 90 90];

b=[aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 ...

aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 ...

aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 ...

aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 ...

aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 ...

aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1];

c=3000*ones(1,108);

ПРИЛОЖЕНИЕ В

(обязательное)

Построение сейсмограмм

Программа в система matlab для построения спектрограмм по трем каналам

%% параметры

% пути к файлам

%p='D:\_____\';% для ноута 11''

p='s:\ATR\';%дляПК

filenameX='2014206055635320_POLE__1_2.atr';

filenameY='2014206055635320_POLE__1_2.atr';

% время

filetime=datenum(2014,07,25,05,56,35.320);

EQtime=datenum(2014,07,25,06,07,02.125);

dt=4;

% фильтр

fmin=3;

fmax=20;

% прочее

delimeter=' ';

nStrok=9;

L=86400;

fs=200;

%% загружаем

% 1-й

data666 = importdata([p filenameX], delimeter, nStrok);

dataX=data666.data;

% 2-й

Продолжение приложения В

data666 = importdata([p filenameY], delimeter, nStrok);

dataY=data666.data;

clear data666;

%% фильтруем

[b,a]=butter(3,[fminfmax]/fs*2);

dataXf=filter(b,a,dataX-mean(dataX));

dataYf=filter(b,a,dataY-mean(dataY));

%% вырезаем

% EQtime - времяземлетрясения

% filetime - времяначалафайла

% они даны в сутках, например 735812.1666240741

% L - кол-во секунд в сутках, =86400

% fs - частота дискретизации, Гц , =200

% dt - длина интересующей нас записи, с

dataXfC=dataXf(((EQtime-filetime)*L*fs-dt*fs/2):((EQtime-filetime)*L*fs+dt*fs/2));

dataYfC=dataYf(((EQtime-filetime)*L*fs-dt*fs/2):((EQtime-filetime)*L*fs+dt*fs/2));

%% ищемразброс

r=max((dataXfC.^2+dataYfC.^2).^.5);

%% вращаем, рисуем

am=99999999999999999999999;

amax=0;

figure(1);

fori=1:18

alf=(i-1)*10;

[ X1,Y1 ] = Func_rotate(dataXfC,dataYfC,alf );

Продолжение приложения В

subplot('position',[.05 .05*(i+.5) .2 .049]);

plot(1:length(X1),X1)

grid on, axis tight, ylim([-r r]);

ylabel(num2str(alf));

ifi>1, set(gca,'xticklabel',''), end

subplot('position',[.3 .05*(i+.5) .2 .049]);

plot(1:length(Y1),Y1)

grid on, axis tight, ylim([-r r]);

ylabel(num2str(alf));

ifi>1, set(gca,'xticklabel',''), end

subplot('position',[.55 .05*(i+.5) .2 .049]);

plot(1:length(X1),X1.^3)

grid on, axis tight, ylim([-r.^3 r.^3]);

ylabel(num2str(alf));

ifi>1, set(gca,'xticklabel',''), end

subplot('position',[.8 .05*(i+.5) .2 .049]);

plot(1:length(Y1),Y1.^3)

grid on, axis tight, ylim([-r.^3 r.^3]);

ylabel(num2str(alf));

ifi>1, set(gca,'xticklabel',''), end

ifstd(Y1)<am

am=std(Y1);

ii=alf;

end

ifstd(Y1)>amax

Продолжение приложения В

amax=std(Y1);

iii=alf;

end

end

set(gcf,'units','normalized','position',[0.15 0.15 .8 .75]);

print(gcf,'-dpng',[p filenameX '.png'],'-r300')

%close(gcf)

% prompt = {'alfa'};

% dlg_title = 'Чемуравноalfa?';

% num_lines = 1;

% def = {num2str(iii)};

% answer = inputdlg(prompt,dlg_title,num_lines,def);

%% вращаем, рисуем ещё раз

am=99999999999999999999999;

amax=0;

figure(2);

fori=1:18

alf=iii-18+(i-1)*2;

[ X1,Y1 ] = Func_rotate(dataXfC,dataYfC,alf );

subplot('position',[.05 .05*(i+.5) .2 .049]);

plot(1:length(X1),X1)

grid on, axis tight, ylim([-r r]);

ylabel(num2str(alf));

ifi>1, set(gca,'xticklabel',''), end

subplot('position',[.3 .05*(i+.5) .2 .049]);

plot(1:length(Y1),Y1)

Продолжение приложения В

gridon, axistight, ylim([-rr]);

ylabel(num2str(alf));

ifi>1, set(gca,'xticklabel',''), end

subplot('position',[.55 .05*(i+.5) .2 .049]);

plot(1:length(X1),X1.^3)

grid on, axis tight, ylim([-r.^3 r.^3]);

ylabel(num2str(alf));

ifi>1, set(gca,'xticklabel',''), end

subplot('position',[.8 .05*(i+.5) .2 .049]);

plot(1:length(Y1),Y1.^3)

grid on, axis tight, ylim([-r.^3 r.^3]);

ylabel(num2str(alf));

ifi>1, set(gca,'xticklabel',''), end

ifstd(Y1)<am

am=std(Y1);

ii=alf;

end

ifstd(Y1)>amax

amax=std(Y1);

iii=alf;

end

end

set(gcf,'units','normalized','position',[0.15 0.15 .8 .75]);

msgbox(num2str(iii))

print(gcf,'-dpng',[p filenameX '__.png'],'-r300')

%close(gcf)

Продолжение приложения В

function [ X1,Y1 ] = Func_rotate( X,Y,alfa )

a=alfa*pi/180;

X1=X*cos(a)-Y*sin(a);

Y1=X*sin(a)+Y*cos(a);

end

ПРИЛОЖЕНИЕ Г

(обязательное)

Расчет азимута землетрясений

Программа для расчета азимута землетрясений, разработанная в системе matlab.

function sp01(fname,fs,wl,overlap, meaning,zmin,zmax)

%load file

delimeter=' ';

nStrok=9;

data666 = importdata(fname, delimeter, nStrok);

data1=data666.data;

clear data666;

% computing

if meaning>1

data1=AntiTrendFast(data1,meaning);

end

[~,F,T,P]=spectrogram(data1,wl,overlap,wl,fs);

F(1:3,:)=[];

P(1:3,:)=[];

%plot 1st

figure;

subplot('position',[0.04 0.75 0.94 0.22]);

set(gca,'fontSize',9)

plot((1:length(data1))./fs,data1); axis tight; grid on;

aaa=get(gca,'ylim');

aaa(1)=aaa(1)-0.02*(aaa(2)-aaa(1));

aaa(2)=aaa(2)+0.02*(aaa(2)-aaa(1));

Продолжение приложения Г

Ylim(aaa);

%plot 2nd

subplot('position',[0.04 0.05 0.94 0.62]);

set(gca,'fontSize',9)

surf(T,F,10*log10(P),'edgecolor','none');

aa1=0; aa2=length(data1)/fs; aa3=F(1,1); aa22=F(size(F)); aa4=aa22(1);

axis ([aa1 aa2 aa3 aa4]);

set(gca,'yscale','log');

colorbar('east');

colormap(jet(4096));

set(gca,'clim',[zminzmax]);

holdon;

a=[.0001 .0001 .0002 .0002 .0003 .0003 .0004 .0004 .0005 .0005 .0006 .0006 .0007 .0007 .0008 .0008 .0009 .0009 ...

.001 .001 .002 .002 .003 .003 .004 .004 .005 .005 .006 .006 .007 .007 .008 .008 .009 .009 ...

.01 .01 .02 .02 .03 .03 .04 .04 .05 .05 .06 .06 .07 .07 .08 .08 .09 .09 ...%9*2*6=108

.1 .1 .2 .2 .3 .3 .4 .4 .5 .5 .6 .6 .7 .7 .8 .8 .9 .9 ...

1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 9 ...

10 10 20 20 30 30 40 40 50 50 60 60 70 70 80 80 90 90];

b=[aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 ...

aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 ...

Продолжение приложения Г

aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 ...

aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 ...

aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 ...

aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1 aa1-1 aa2+1 aa2+1 aa1-1];

c=3000*ones(1,108);

plot3(b,a,c,':b'); hold off;