%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% do_averages.m
% This program calculates mean intensitites from data extracted from the
% different logs (mostly painful hand-extraction from totally incompatible 
% and non-standard EXCEL files). This program takes A LONG time to run as
% 
% This program saves 3 files with data that will be useful for plotting and
% other analyses:
%
%   U1309D_averages5m   5-m averages
%   U1309D_arrays1cm    1 cm arrays for the different 'parameters'
%   U1309D_cleandata    Cleaned, continuous U1309D intervals
%
%   Files contain veining, cataclasis, open fractures, drilling,
%   serpentinite foliation, magmatic and plastic intensities + lithology
%
%   Calculations use:
%       downhole.m - to do averages of a given depth interval
%       holearray.m - creates arrayfrom intervals (i.e., at 1 cm spacing)
%       dodepth.m - gives depths (shifted & curated) from core, section and
%       depth from top of section
%
% Javier Escartin - 27/2/05 - escartin@ipgp.jussieu.fr
%
% If files are updated, you'll need to re-run this program and modify all
% subsequent pltos and calculations.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;clf; warning off
t=cputime

display('loading data')

% Loading brittle log from Expedition 304 and 305
load U1309D_brittlelog.dat;brittle305D=U1309D_brittlelog;
load U1309Dtop_brittlelog.dat; brittle304D=U1309Dtop_brittlelog;
clear U1309D_brittlelog U1309Dtop_brittlelog

% Separate the different intensity arrays for more 'comfortable' analysis
doi=[brittle304D(:,[1 2 3]);brittle305D(:,[1 2 3])]; % Open cracks 0-3
ddi=[brittle304D(:,[1 2 4]);brittle305D(:,[1 2 4])]; % drill fractures 0-3
dsi=[brittle304D(:,[1 2 5]);brittle305D(:,[1 2 5])]; % Serpentinization foliation  0-3
dci=[brittle304D(:,[1 2 6]);brittle305D(:,[1 2 6])]; % cataclastic  0-5

% Vein - these have been consolidated on one of the Structure tables in
% 305, while 304 had a dedicated vein table
load U1309Dup_veins.dat;jk=U1309Dup_veins; clear U1309Dup_veins;

% The up veins do not have a curated+shifted depth - so let's do it
load U1309D_shift.dat;s= U1309D_shift; % this is the shift table
load U1309D_cores.dat; d=U1309D_cores; %this is the curated depth for each section

% calculate depth of top and bottom of each interval (uses dodepth.m). For
% Upper part only (tables do not contain curated/shifted depths
nt=dodepth(jk(:,1:3),d,s);nt=nt(:,4); % top of interval
nb=dodepth(jk(:,[1:2 4]),d,s);nb=nb(:,4); % bottom of interval
dvi1=[nt nb jk(:,5)]; % creating a regular depth top depth bottom value array (3 cols)
i=find(dvi1(:,3)>=0);dvi1=dvi1(i,:); % remove -1 flags (heritage from Expedition 304)
i=find(dvi1(:,3)>5);dvi1(i,3)=dvi1(i,3)*5./dvi1(i,3); % remove values larger than 5 (should not exist, but...)
clear jk

% Load another vein datafile extracted from 304 EXCEL tables
load U1309Dmid_veinint.dat
dvi0=U1309Dmid_veinint;clear U1309Dmid_veinint
i=find(dvi0(:,3)>=0);dvi0=dvi0(i,:); % Do some clean up as above
i=find(dvi0(:,3)>5);dvi0(i,3)=dvi0(i,3)*5./dvi0(i,3);

% Join arrays for a single U1309D vein dataset
dvi=[dvi1;brittle305D(:,[1 2 7])];
clear dvi0 dvi1

clear brittle304D brittle305D

% high T data loading
load U1309D_magmaticint.dat;dmi=U1309D_magmaticint;
load ../Data_clean/U1309D_plasticint.dat;dpi=U1309D_plasticint;

% clean for -1 stupid numbers (if not done above) for all arrays
i=find(dvi(:,3)>=0);dvi=dvi(i,:);% veining
i=find(doi(:,3)>=0);doi=doi(i,:);% open cracks
i=find(ddi(:,3)>=0);ddi=ddi(i,:);% drilling
i=find(dsi(:,3)>=0);dsi=dsi(i,:);% serpentinite foliation
i=find(dci(:,3)>=0);dci=dci(i,:);% cataclasis
i=find(dmi(:,3)>=0);dmi=dmi(i,:);% magmatic foliation
i=find(dpi(:,3)>=0);dpi=dpi(i,:);% plastic foliation
% Load the lithological column
!awk '{print $1, $2, $4}' U1309D_lithunit.dat > lithomat.dat
load lithomat.dat
l=lithomat;clear lithomat;litho=l;!\rm lithomat.dat;
save U1309D_cleandata dvi doi ddi dsi dci dmi dpi litho % save the interval data
clear litho

% Do averages - here we define the array - 0-1420 mbsf with a 5 m averaging
% window
zb=0;ze=1420;dz=5;

display('start the calculations of averages')
Dvein=downhole(dvi, zb, ze, dz);display('done veins U1309D')
Dopen=downhole(doi, zb, ze, dz);display('done open fractures U1309D')
Ddril=downhole(ddi, zb, ze, dz);display('done drill fractures 304D')
Dserp=downhole(dsi, zb, ze, dz);display('done serpentinization foliation U1309D')
Dcata=downhole(dci, zb, ze, dz);display('done cataclasis U1309D')
Dmagm=downhole(dmi, zb, ze, dz);display('done magmatic foliation 304D')
Dplas=downhole(dpi, zb, ze, dz);display('done plastic foliation 304D')
save U1309D_averages5m Dvein Dopen Ddril Dserp Dcata Dmagm Dplas % save the average data

% Do arrays for calculations - here we define the array - 0-1420 mbsf with
% a 1 cm array - no overlap of sections
zb=0;ze=1420;dz=.1;
display ('do arrays at 1 cm')
Hvein=holearray(dvi, zb, ze, dz);display('done veins U1309D')
Hopen=holearray(doi, zb, ze, dz);display('done open fractures U1309D')
Hdril=holearray(ddi, zb, ze, dz);display('done drill fractures 304D')
Hserp=holearray(dsi, zb, ze, dz);display('done serpentinization foliation U1309D')
Hcata=holearray(dci, zb, ze, dz);display('done cataclasis U1309D')
Hmagm=holearray(dmi, zb, ze, dz);display('done magmatic foliation 304D')
Hplas=holearray(dpi, zb, ze, dz);display('done plastic foliation 304D')
Hlith=holearray(l  , zb, ze, dz);display('done lithology array')
save U1309D_arrays1cm Hvein Hopen Hdril Hserp Hcata Hmagm Hplas Hlith

cputime-t