clear;clf;

%load data\vars_Harvey.mat;         %22
%load data\vars_Michael.mat;        %34
load data\vars_Patricia.mat;        %7   
%load data\vars_Rita.mat;           %50
%load data/vars_20181010;

E=6.11*exp(2500000/416.*ta./273.16./(ta+273.16));
e=6.11*exp(2500000/416.*td./273.16./(td+273.16));

rh=e./E;

q=0.622*e./(p-e);

th=(ta+273.16).*(1000./p).^(287/1004);
th_e=th.*(1+2500000./(1004*(ta+273.16)).*q);

[d1 d2]=size(th_e);

sb(1:10,1:d2)=nan;
se(1:10,1:d2)=nan;

%d2=22;

for n=1:d2   %1
n
rh1=rh(:,n);

[I]=find(rh1>0.93);
  if isempty(I)  %2 
  continue
  else           %2
    nb=1;
    ne=find_sat(nb,I);
    sb(1,n)=I(nb);
    se(1,n)=I(ne);
    nb=ne+1;
    if nb > length(I)
        continue
    end
    ne=find_sat(nb,I);
    sb(2,n)=I(nb);
    se(2,n)=I(ne);
    
    nb=ne+1;
    if nb > length(I)
       continue 
    end
    ne=find_sat(nb,I);
    sb(3,n)=I(nb);
    se(3,n)=I(ne);
    
    nb=ne+1;
    if nb > length(I)
        continue
    end
    ne=find_sat(nb,I);
    sb(4,n)=I(nb);
    se(4,n)=I(ne);
    
    nb=ne+1;
    if nb > length(I)
        continue
    end
    ne=find_sat(nb,I);
    sb(5,n)=I(nb);
    se(5,n)=I(ne);

    nb=ne+1;
    if nb > length(I)
        continue
    end
    ne=find_sat(nb,I);
    sb(6,n)=I(nb);
    se(6,n)=I(ne);

    nb=ne+1;
    if nb > length(I)
        continue
    end
    ne=find_sat(nb,I);
    sb(7,n)=I(nb);
    se(7,n)=I(ne);

    nb=ne+1;
    if nb > length(I)
        continue
    end
    ne=find_sat(nb,I);
    sb(8,n)=I(nb);
    se(8,n)=I(ne);

    nb=ne+1;
    if nb > length(I)
        continue
    end
    ne=find_sat(nb,I);
    sb(9,n)=I(nb);
    se(9,n)=I(ne);

    nb=ne+1;
    if nb > length(I)
        continue
    end
    ne=find_sat(nb,I);
    sb(10,n)=I(nb);
    se(10,n)=I(ne);

  end           %2

end      %1

k=1;
for n=1:d2 %1
th1=th(:,n);
th_e1=th_e(:,n);
q1=q(:,n);

  var=sb(:,n);
  var1=var(abs(var)>=0);
  var=se(:,n);
  var2=var(abs(var)>=0);
  var3=var2-var1;
  var1=var1(var3>=2 & var1>3);
  var2=var2(var3>=2 & var2<297);
  for i=1:length(var1) %2
  dth_e(k)=nanmean(th_e1(var1(i)-2:var1(i)-1))-nanmean(th_e1(var1(i):var1(i)+1));
  dq(k)=(nanmean(q1(var1(i)-2:var1(i)-1))-nanmean(q1(var1(i):var1(i)+1)))*mean(th1(var1(i)-2:var1(i)+1))*2;
  k=k+1;
  end %2
  for i=1:length(var2) %2
  dth_e(k)=nanmean(th_e1(var2(i)+1:var2(i)+2))-nanmean(th_e1(var2(i)-1:var2(i)));
  dq(k)=(nanmean(q1(var2(i)+1:var2(i)+2))-nanmean(q1(var2(i)-1:var2(i)+1)))*mean(th1(var2(i)-1:var2(i)+2))*2;
  k=k+1;
  end %2
end %1

%save data/Rita_ent_dat.mat dth_e dq;
%save data/Harvey_ent_dat.mat dth_e dq;
%save data/Michael_ent_dat.mat dth_e dq;
save data/Patricia_ent_dat.mat dth_e dq;

plot(dth_e,dq,'*');
axis([-12,0,-12,0]);
hold on;
plot([-12 0],[-12 0],'k');
hold off;

