% EPS 131 : HW-02, question 1 
% Laure Zanna
% ---------------------
 clear all; close all;
 
 % ------------------------------------------------------------------
 % STEP 1: select data (in binary random access form) from the iri website at
 % (http://iridl.ldeo.columbia.edu/SOURCES/.LEVITUS94/.ANNUAL/.temp/)   
 % I called the 2 files data_eq.r4 and data_45.r4
 
 % ------------------------------------------------------------------
 % STEP 2: Use the f77read.m file in order to read the binary data files
 % read the data  for the temperature [X,Y=longitude,Latitude Z=depth] 
 % Z=[33 levels] 
 % --------------------------------------------------------------------
 % Equatorial Pacific at Lat=0.5S, Long=179.5E
 T_eq=f77read('data_eq.r4',[1 33],'real*4','b');
 % Pacific at Lat=45.5N, Long=160.5W
 T_45=f77read('data_45.r4',[1 33],'real*4','b');
 % define the depth
 depth=[0.0 10 20 30 50 75 100 125 150 200 250 300 400 500 600 700 800 900 ...
        1000 1100 1200 1300 1400 1500 1750 2000 2500 3000 3500 4000 4500 5000 5500];
 
 % --------------------------------------------------------------------
 % STEP 3: 
 % plot the 2 verical profiles as function of depth
 % --------------------------------------------------------------------
 %define x and y labels for plots
 labels=char('Depth [m]','Temperature [\circ C]');
 
 figure(1)
 subplot(2,1,1); plot(depth,T_eq,'.-')
 title('Equatorial Pacific Lat=0.5S , Long=179.5E','fontsize',14)
 ylabel(labels(2,:),'fontsize',13)
 
 subplot(2,1,2) ; plot(depth,T_45,'.-')
 title('Pacific at Lat=45.5N , Long=160.5W','fontsize',14)
 xlabel(labels(1,:),'fontsize',13); ylabel(labels(2,:),'fontsize',13)

 
 % --------------------------------------------------------------------
 % STEP 4: fitting the data from 500m to the bottom of the ocean
 % --------------------------------------------------------------------
 % (1) plot from 500m to the bottom of the ocean
 
 figure(2)
 subplot(1,2,1); plot(depth(14:33),T_eq(14:33),'.-');hold on;
 title('Equatorial Pacific from 500m','fontsize',14)
 xlabel(labels(1,:),'fontsize',13); ylabel(labels(2,:),'fontsize',13); 
 
 subplot(1,2,2) ; plot(depth(14:33),T_45(14:33),'.-');hold on;
 title('Pacific at Lat=45.5N from 500m','fontsize',14)
 xlabel(labels(1,:),'fontsize',13); ylabel(labels(2,:),'fontsize',13);
 
 
 % (2) The data cannot contain any NaN in order to fit, therefore 
 % we remove the last point at depth of 5500m (its a NaN),  
 
 % select the data from 500m to 5000m 
 depth=depth(14:32)'; T_eq=T_eq(14:32); T_45=T_45(14:32);
 
 % (3) define the baseline from the plots : the baseline in our case is the nearly 
 % constant temperature reached as we get closer to the bottom of the ocean
 % We call them T_0_eq and T_0_45
  T_0_eq=1.2585;
  T_0_45=1.517;
 
 % (4) We wish to fit our profiles to an exponential such that
 % T=T_0+a*exp(-z/H)  ->   T-T_0= a*exp(-z/H)
 % We can transorm this function to a polynomial by taking the log:
 % ln(T-T_0) = ln (a*exp(-z/H))= ln(a) + (-z/H)
 % We can now fit ln(T-T_0) = ln(a) + (-z/H) 
 
 %fitting to polynomial of order 1 : f1(x) = p1*x + p2
 f1_eq = fit(depth,log(T_eq-T_0_eq)','poly1')
 f2 = fit(depth,log(T_45-T_0_45)','poly1')
 
 % (5) : use the results from Matlab fitting the recontruction 
 % the exponential curve
 % we have ln(a)=p1 -> a=exp(p2) and -1/H=p1;
 
 %at the equator
 p1_eq =-0.00213; p2_eq =3.542 ;
 H_eq=-1.5*1/p1_eq; a_eq=0.4*exp(p2_eq);
 % at 45N
 p1_45 = -0.001343; p2_45 =1.682;
 H_45=-1/p1_45; a_45=exp(p2_45); 
 
 % (6) Plot the exponential curves on top of the data...
 fit_eq=T_0_eq+a_eq*exp(-depth/H_eq);
 figure(2); hold on;
 subplot(1,2,1); plot(depth,fit_eq,'g-','linewidth',2);legend('data','LS fit')
 fit_45=T_0_45+a_45*exp(-depth/H_45);
 figure(2); hold on;
 subplot(1,2,2); plot(depth,fit_45,'g-','linewidth',2);legend('data','LS fit')
