%% Master script pure equations
% Written May 2007 (kieran Helm)
clear; clc; 

d=[230:2:258 260:284];
dep=26;

%% Load data
[t1,t2,w_lat,w_lon]=kh_load_nddata('ns4',d,'t'); %t1=historical, t2=modern
[p1,p2,w_lat,w_lon]=kh_load_nddata('ns4',d,'p');
[s1,s2,w_lat,w_lon]=kh_load_nddata('ns4',d,'s');
comment=['loaded data ']

%% Remove anything <100m
% a=find(p1<100); p1(a)=NaN; s1(a)=NaN; t1(a)=NaN; 
% a=find(p2<100); p2(a)=NaN; s2(a)=NaN; t2(a)=NaN; 

%% Get subduction coeff (use in situ temp)
% [Aw Af Ah Aw_r2 Af_r2 Ah_r2]=kh_subduction_coeff(s1,s2,t1,t2,p1,p2); 
% save('subduction_coeff','Aw','Af','Ah','Aw_r2','Af_r2','Ah_r2');
load('subduction_coeff')

%% Get terms
[tz tn sz sn Ntz Nsz rp]=kh_get_all_subduction_terms(s1,s2,t1,t2,p1,p2);

%% Plot subduction proportions
% kh_plot_subduction_proportions(Aw_r2,Af_r2,Ah_r2,w_lat,w_lon,d,37,'squared')

%% Other plots
% a=find(w_lat>-90);
% a=find(w_lat>-40 & w_lat<-25 & w_lon>30 & w_lon<115);
% gAw=kh_interp_ignore(Aw_r2(:,a),w_lat(a),w_lon(a),[ones(40*1)*100],5,[5 10],'3D');
% gAf=kh_interp_ignore(Af_r2(:,a),w_lat(a),w_lon(a),[ones(40*1)*100],5,[5 10],'3D');
% gAh=kh_interp_ignore(Ah_r2(:,a),w_lat(a),w_lon(a),[ones(40*1)*100],5,[5 10],'3D');
% plot(d,kh_nanmean2(gAw,3),'r-*'); hold on
% plot(d,kh_nanmean2(gAf,3),'b-*'); hold on
% plot(d,kh_nanmean2(gAh,3),'g-*'); hold on; grid on