% 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(1,2,1); plot(depth,T_eq,'.-');hold on;
 title('Equatorial Pacific Lat=0.5S , Long=179.5E','fontsize',14)
 xlabel(labels(1,:),'fontsize',13); ylabel(labels(2,:),'fontsize',13)
 
 subplot(1,2,2) ; plot(depth,T_45,'.-');hold on;
 title('Pacific at Lat=45.5N , Long=160.5W','fontsize',14)
 xlabel(labels(1,:),'fontsize',13); ylabel(labels(2,:),'fontsize',13)

 T_0_eq=1.2585+0*depth;
  T_0_45=1.5170+0*depth;
  
  figure(1)
  subplot(1,2,1); plot(depth,T_0_eq,'-');
  subplot(1,2,2) ; plot(depth,T_0_45,'-');hold on;
 
 % --------------------------------------------------------------------
 % 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+0*depth;
  T_0_45=1.5170+0*depth;
  
  figure(1)
  subplot(1,2,1); plot(depth,T_0_eq,'-');
  subplot(1,2,2) ; plot(depth,T_0_45,'-');hold on;
  
  % use 2 other point on the graph
  % at depth=500m, T_eq=8.2359; T_45= 4.2564
  % T=T_0+a*exp(-z/H)
  % eq: 8.2359=1.2585+ a*exp(-500/H) -> 6.9774=a*exp(-500/H)
  % 45N: 4.2564=1.5170+ a*exp(-500/H) -> 2.7394=a*exp(-500/H)
  
  % z=500+H, at the eq: T(500+H)=3.8253C; at 45N T(500+H)=2.5248
  
  a_eq=14.2529; a_45=5.1179;
  H_eq=700 ; H_45=800;
  T_0_eq=1.2585; T_0_45=1.5170;
  % (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','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','fit')

