% This program estimates the density of the rocks with velocity as input
% (example: Baikal Rift)

% Estimate the depth of each crustal layer (Christensen and Mooney, 1995)

DepthUC=[15,18,19];
DepthMC=[28,31,30];
DepthLC=[41,45,49];

% Define the depth of sediments and mantle
DepthSed=zeros(1,3);
DepthMantle(1:3)=150;

DepthCrust=[DepthSed;DepthUC;DepthMC;DepthLC;DepthMantle];

% Define the velocity of each crustal layer

Vuc1=[6.0,6.3;6.1,6.4;6.1,6.4];
Vmc1=[6.4,6.7;6.6,6.8;6.5,6.7];
Vlc1=[7.2,7.3;7.4,7.6;7.2,7.3];

% Estimate the velocity of each layer

Vuc=mean(Vuc1');
Vmc=mean(Vmc1');
Vlc=mean(Vlc1');

% Define the coefficients of the velocity-density relation

a=[4929;5055;5141;5212;5281];
b=[-13294;-14094;-14539;-14863;-15174];
DEPTH=[10:10:50];
DepthRocks=[DepthUC./2;((DepthMC-DepthUC)./2)+DepthUC;((DepthLC-DepthMC)./2)+DepthMC];

% Define the coefficients by linear interpolation for each depth

% upper crust

auc=interp1(DEPTH,a,DepthRocks(1,:),'linear','extrap');
buc=interp1(DEPTH,b,DepthRocks(1,:),'linear','extrap');

% middle crust

amc=interp1(DEPTH,a,DepthRocks(2,:),'linear','extrap');
bmc=interp1(DEPTH,b,DepthRocks(2,:),'linear','extrap');

% lower crust

alc=interp1(DEPTH,a,DepthRocks(3,:),'linear','extrap');
blc=interp1(DEPTH,b,DepthRocks(3,:),'linear','extrap');

% Estimate the density of each crustal layer

% Upper crust

Rouc1=auc+buc./Vuc;

% Middle crust

Romc1=amc+bmc./Vmc;

% Lower crust

Rolc1=alc+blc./Vlc;

% Define the density of sediments and mantle
DenSed=zeros(1,3);
DenM(1:3)=3350;

RoRocks=[DenSed;Rouc1;Romc1;Rolc1;DenM];

save DensityRocks.dat RoRocks -ASCII
save DepthCrust.dat DepthCrust -ASCII
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Estimate the density of sediments (Brocher, 2005)

VpSed=[1.9,3.6;4.5,5.3];
Vp1=mean(VpSed');
Vp1=mean(Vp1);

DenSed1=1.6612*Vp1-0.4721.*Vp1.^2 + 0.0671.*Vp1.^3-0.0043.*Vp1.^4+0.000106.*Vp1.^5;