disp('timatch_3_1')
fclose('all');

skip_it=0;
if skip_it==0
clear all;
jaz=zeros(16,1440);
jazx=zeros(6,1440);
jazy=zeros(13,1440);

jd=input('julian day, ddd   ','s');

f=['D:\hp_ascii\proc_son' jd  '.txt'];
flist=fopen(f,'r');
jazx=fscanf(flist,'%g%g%g%g%g%g',[6,inf]);

e=['D:\hp_ascii\proc_scs' jd  '.txt'];
flist2=fopen(e,'r');
jaz=fscanf(flist2,'%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g',[16,inf]);

g=['D:\hp_ascii\proc_hp' jd  '.txt'];
flist3=fopen(g,'r');
jazy=fscanf(flist3,'%g%g%g%g%g%g%g%g%g%g%g%g%g',[13,inf]);

end;

tm_in=[124.00 129.10 130.083 131.14 132.10];			%times for snake in and out of water
tm_ot=[128.60 129.50 130.158 131.27 158.00];
zsnk=ones(1,288);ztsg=ones(1,288);

tm_s=[124.00    149.71133];								%time entering S lat.
tm_n=[126.4925  158.00];

[ax,bx]=size(jazx);												%jazx
[c,d]=size(jaz);												%jaz
[ay,by]=size(jazy);
tn=60*24;														%number of minutes in a day
n=1:tn;
jad0=str2num(jd)+(n-1)/24/60;jad0(tn)=jad0(tn)+1e-5;		%	don't ask
jadson=jazx(1,1:bx);
jadscs=jaz(1,1:d);
jadmns=jazy(1,1:by);

for ii=2:length(jadscs);
   if jadscs(ii)<0,jadscs(ii)=jadscs(ii-1)+.000001;end;
end;

%%%%%%%%**********   sonic loop
clear jazx1;
%fill in beginning empty spots 
jk=1;
while jad0(jk)<jadson(1),
   jazx1(:,jk)=jazx(:,1);
   jk=jk+1;
end;

for i=jk:tn;
   t=find(jadson>jad0(i));j=t(1);
   dt=jadson(j)-jadson(j-1);Dt=jad0(i)-jadson(j-1);
   jazx1(:,i)=jazx(:,j-1)*(1-Dt/dt)+jazx(:,j)*Dt/dt;
   if j==bx;break;end;
 end;
jk=tn;
while jad0(jk)>=jadson(bx),
   jazx1(:,jk)=jazx(:,bx);
   jk=jk-1;
end;

%%%%%%%%*********  means loop
clear jazy1;
%fill in beginning empty spots 
jk=1;
while jad0(jk)<jadmns(1),
   jazy1(:,jk)=jazy(:,1);
   jk=jk+1;
end;

for i=jk:tn;
   t=find(jadmns>jad0(i));
   if length(t)>0,
      j=t(1);
   else
      j=2;
   end;
   dt=jadmns(j)-jadmns(j-1);Dt=jad0(i)-jadmns(j-1);
   jazy1(:,i)=jazy(:,j-1)*(1-Dt/dt)+jazy(:,j)*Dt/dt;
   if j==by;break;end;
end;
jk=tn;
while jad0(jk)>=jadmns(by),
   jazy1(:,jk)=jazy(:,by);
   jk=jk-1;
end;

%%%%%%%%%%%%%****************   ship data loop
%first setup to interpolate fricking angles
cog=jaz(5,:);
sog=jaz(6,:);
lrg=jaz(13,:);
imu=jaz(15,:);
imd=jaz(16,:);

rdcon=pi/180;
jaz(5,:)=sog.*cos(cog*rdcon);			%sogx
jaz(6,:)=sog.*sin(cog*rdcon);			%sogy
jaz(15,:)=imu.*cos(imd*rdcon);
jaz(16,:)=imu.*sin(imd*rdcon);
lrgc=cos(lrg*rdcon);
lrgs=sin(lrg*rdcon);

clear jaz1;
												%fill in beginning empty spots 
jk=1;
while jad0(jk)<jadscs(1),
   jaz1(:,jk)=jaz(:,1);
   lrgcm(jk)=lrgc(1);
   lrgsm(jk)=lrgs(1);
   jk=jk+1;
end;

for i=jk:tn;
   t=find(jadscs>jad0(i));
   if~isempty(t)
	   j=t(1);
   	dt=jadscs(j)-jadscs(j-1);Dt=jad0(i)-jadscs(j-1);
   	jaz1(:,i)=jaz(:,j-1)*(1-Dt/dt)+jaz(:,j)*Dt/dt;
   	lrgcm(i)=lrgc(j-1)*(1-Dt/dt)+lrgc(j)*Dt/dt;
   	lrgsm(i)=lrgs(j-1)*(1-Dt/dt)+lrgs(j)*Dt/dt;
	else
   	break;
	end;
end;
jk=tn;
while jad0(jk)>=jadscs(d),
   jaz1(:,jk)=jaz(:,d);
   lrgcm(jk)=lrgc(d);
   lrgsm(jk)=lrgs(d);
   jk=jk-1;
end;

%************************   do 30 minute averages
jax=1440;
jk=288;
ode=jaz1(14,:);
sodx=ode.*lrgcm;
sody=ode.*lrgsm;
uu=jazx1(2,:)/.95;
vv=jazx1(3,:)/1.15;
un=uu.*lrgcm+vv.*lrgsm;
ue=uu.*lrgsm-vv.*lrgcm;

for j=1:jk
  	jl=5*(j-1)+1;jj=min(jax,5*j);
     jazx5(:,j)=mean(jazx1(:,jl:jj)')';
     jazy5(:,j)=mean(jazy1(:,jl:jj)')';
     jaz5(:,j)=mean(jaz1(:,jl:jj)')';
     lrgc5(j)=mean(lrgcm(jl:jj)')';
     lrgs5(j)=mean(lrgsm(jl:jj)')';
     sodx5(j)=mean(sodx(jl:jj)')';
     sody5(j)=mean(sody(jl:jj)')';
     un5(j)=mean(un(jl:jj)')';
     ue5(j)=mean(ue(jl:jj)')';
     vrsp(j)=var(jaz1(5,jl:jj))+var(jaz1(6,jl:jj));
     jad5(j)=jad0(jl);
end;
   
%*********************  then uncomponentize angles

U=jazx5(2,:);										%.95 is empirical distortion factor
V=jazx5(3,:);										%1.15 is other empirical distortin factor
w=jazx5(4,:);
reldir=mod(atan2(V,-U)/rdcon+360,360);		%signs convert jeff coord to conventional
relsp=sqrt(U.*U+V.*V);
sogx=jaz5(5,:);
sogy=jaz5(6,:);
imux=jaz5(15,:);
imuy=jaz5(16,:);

sogm5=sqrt(sogx.*sogx+sogy.*sogy);					%1-min sog
cogm5=mod(atan2(sogy,sogx)/rdcon+360,360);
imum5=sqrt(imux.*imux+imuy.*imuy);					%1-min imet speed, m/s
imdm5=mod(atan2(imuy,imux)/rdcon+360,360);
lrgm5=mod(atan2(lrgs5,lrgc5)/rdcon+360,360);			%laser ring gyro
odec5=sqrt(lrgc5.*lrgc5+lrgs5.*lrgs5);

%**********  do some true wind calcs
un1=un5+sogx;
ue1=ue5+sogy;
un2=un5+sodx5;
ue2=ue5+sody5;
s1=sqrt(un1.*un1+ue1.*ue1);
s2=sqrt(un2.*un2+ue2.*ue2);
s3=sqrt(un1.*un1+ue1.*ue1+w.*w);

dir1=mod(atan2(-ue1,-un1)*180/pi+360,360);
dir2=mod(atan2(-ue2,-un2)*180/pi+360,360);

%****************   sea temps
ts_snk=jazy5(10,:);
ts_tsg=jaz5(7,:);
sal_tsg=jaz5(8,:);
ts=ts_tsg;zt=ztsg*5;
[n,m]=size(tm_in);
jg=[];
for k=1:m
   jf=find(jad5>tm_in(k) & jad5<tm_ot(k));%find when snake is in water
   jg=[jg jf];
end;
ts(jg)=ts_snk(jg);							%water temperature if snake is in
zt(jg)=zsnk(jg)*.05;							%depth of water temp measurement
qs=qsea(ts);

%%%%%************   air temps
ta_etl=jazy5(11,:);
ta_im=jaz5(10,:);
ta=min(ta_etl,ta_im);

%***************    relative hums
rh_etl=jazy5(12,:);
rh_im=jaz5(9,:);
qa_etl=qair([ta_etl' rh_etl']);
qa_im=qair([ta_im' rh_im']);
x=[qa_im qa_etl];
qa=qa_etl';												%qa=mean(x');
%**************   solar fluxes
psp1=jazy5(4,:);
psp2=jazy5(5,:);
psp_im=jaz5(11,:);
x=[psp1' psp2'];
psp=mean(x');

%**************  IR fluxes
pir1=jazy5(2,:);
pir2=jazy5(3,:);
x=[pir1' pir2'];
pir=mean(x');

%**************************
org=max(0,jazy5(13,:)-.3);
lat=jaz5(3,:);
lon=jaz5(4,:);
[n,m]=size(tm_s);
jg=[];
for k=1:m
   jf=find(jad5>tm_s(k) & jad5<tm_n(k));					%find when south lat
   jg=[jg jf];
end;
lat(jg)=-lat(jg);

shp_spd=odec5;
shp_hed=lrgm5;

sig_sp=sqrt(vrsp);


prt_f=0;
prt_f=input('append efb data file? 0=n, 1=y   ');
np=jk;
if prt_f==1;
  prt_jas_efb;
end;

for k=1:jk
   %   u    0  ts  ta      qs     qa   rs      rl   rain     zi  p   zu    zt   zq
   x=[s1(k) 0 ts(k) ta(k) qs(k) qa(k) psp(k) pir(k) org(k) 600 1010 17.5 15.5 15.5];

   y=cor26(x);
   hsb(k)=y(1);
   hlb(k)=y(2);
   taub(k)=y(3);
   tcol(k)=y(11);
   qcol(k)=y(12);
   rf(k)=y(14);
   usr(k)=y(8);
   zo(k)=y(4);
end;
Rns=psp1*.955;
Rnl=.97*(pir1-5.67e-8*(ts+273.15).^4);
hnet=Rns+Rnl-hsb-hlb-rf;
hm=mean(hnet);sm=mean(hsb);lm=mean(hlb);Sm=mean(Rns);Lm=mean(Rnl);
figure;plot(jad5,hnet);xlabel('Julian Day (1999)');ylabel('Net Heat Flux (W/m^2)');Title(['Heat Fluxes:  Net= ' num2str(hm) ' Hl= ' num2str(lm)  ' Hs= ' num2str(sm) ' Rs= ' num2str(Sm) ' Rl= ' num2str(Lm)]);
figure;plot(jad5,s1,jad5,s2);xlabel('Julian Day (1999)');ylabel('True Wind Speed (m/s)');Title('Rel. Earth=blu; Rel. Ocean=grn');
figure;plot(jad5,hsb, jad5, hlb,jad5, Rnl,jad5,rf,'.');xlabel('Julian Day (1999)');ylabel('Heat Flux Components (W/m^2)');Title('Hs=blu; Hl= grn; Rnl=red; Rf=dots');

prt_f=0;
prt_f=input('append flux data file? 0=n, 1=y   ');
np=jk;
if prt_f==1;
  prt_jas_flux;
end;
