disp('read_shp_AMMA_2008.m')
fclose('all');

if ddd>99,
   jd=num2str(ddd);
else   
   jd=['0',num2str(ddd)];
end;

g=0.0098;							%adiabatic lapse rate
eps_w=0.97;							%emissivity of water
sig_sb=5.67e-8;					%Stefan Boltzmann constant
blk=1;								%10-min index 
jax=1;
jaxx=1;
backchk=0;
clear st1 st2 st3;
jam=hhh;

if jam<10,
   hr=['0',num2str(jam)];
else
   hr=num2str(jam);
end;

dfl=['C:\Data\AMMA\day' ,jd,'\scs008' jd hr '_raw.txt'];
disp(dfl)

st2=[];temp=[];
flist=fopen(dfl,'r');
%  dmp2=fgets(flist);
if flist>0,   					%if the file exists, 
   while feof(flist)==0,

    try
     temp =textscan(flist,'%2f%2f%3f %f %f%c %f%c %f %f %f %f %f %f %f %f %f %f %f %f %f %*[^\n]','delimiter', ',', 'headerlines', 1,  'emptyvalue', NaN);
     temp{1,6}=double(temp{1,6});
     temp{1,8}=double(temp{1,8});
     st2=[st2 cell2mat(temp)'];
     
    catch 
%      for ll=2:21
%          if length(temp{1,ll})~=length(temp{1,1})
%              temp{1,ll}(length(temp{1,1}))=NaN;
%          end;
%      end;

        for ll=1:21
             [l21,c21]=size(temp{1,21});
             if length(temp{1,ll})~=l21
                 temp{1,ll}(l21+1)=[];
             end;
         end;
        if isempty(temp{1,1});
            temp{1,9}=temp{1,9}';
        end;
     st2=[st2 cell2mat(temp)'];
    end;    
   
   end;  %end while
   
if ddd<75
    st2(18,:)=0; %no odec data
% elseif ddd>=75
%     st2(:,any(isnan(st2),1))=[]; %remove lines containing NaNs (bad lines, too long...)
end;    

    st2(:,any(isnan(st2),1))=[]; %remove lines containing NaNs (bad lines, too long...)

    
else
     st1=ones(21,1800).*NaN;
     st2=ones(21,1800).*NaN;
end; %end if

rdcon=pi/180;

tgyroref=(0:2:3599);
tgyro=st2(1,:)*60+st2(2,:)+st2(3,:)/1000; %real raw time in s
% tmotion=(0:1/20:71999/20);
tmotion=(0:1/10:35999/10);

%10/23/06 - Correction of the ship sampling error. 
%Before the routine was correcting the cog, sog and lrg in order to have samples of 1800 values(on line 178 and after). 
%The problem was that the original data were 1682 values. Thus
%an error was introduced at the end...
%size(st2) = [17,1682];
%size(st2t) = [1800,17];    AS if we had 30 data/min... 

% tt=0:60/(length(st2)/60):3599;     %true sampling time from the ship (=2.1403 s)
st2t=interpft(st2',1800);       %interpolation by FFT on a 2s sampling time.
% st2t=interp1(tgyro,st2',tgyroref,'nearest');       %interpolation to a 2s sampling time.
st2t=st2t';

% figure;plot(tgyro,st2(9,:),'.',tgyroref,st2t(9,:),'g.')

%replace NaN values by 1st good real points
% ii=find(isnan(st2t));
% if ~isempty(ii)
%     if ddd>=75
%         for jj=1:length(ii)/21
%             st2t(:,jj)=st2t(:,(length(ii)/21)+1);
%         end;
%     elseif ddd<75
%         for jj=1:length(ii)/20
%             st2t(:,jj)=st2t(:,(length(ii)/20)+1);
%         end;
%     end;
% end;

% if ii(length(ii))<=20
%     st2t(:,1)=st2t(:,2);
% end;

jax2=length(st2t);

cog=st2t(9,1:jax2);     %in degrees
sog=st2t(10,1:jax2);     %in knots

lrg=st2t(17,1:jax2);    %in degrees

% if ddd>=75
    odec =    st2t(18,:);
    imu =     st2t(19,:);   
    imd =     st2t(20,:);
    ir  =     st2t(21,:);
% elseif ddd<75
%     odec =    ones(1,length(st2t))*NaN;
%     imu =     st2t(18,:);   
%     imd =     st2t(19,:);
%     ir  =     st2t(20,:);
% end;
                
% imu=st2t(16,1:jax2);
% imd=st2t(17,1:jax2);
jk=ceil(jax2/30);
% st3=zeros(18,jk);
% rdcon=pi/180;
st2t(9,1:jax2)=sog.*cos(cog*rdcon);
st2t(10,1:jax2)=sog.*sin(cog*rdcon);

% if ddd>=75
    st2t(19,1:jax2)=imu.*cos(imd*rdcon);
    st2t(20,1:jax2)=imu.*sin(imd*rdcon);
    st3=zeros(21,jk);
% elseif ddd<75
%     st2t(18,1:jax2)=imu.*cos(imd*rdcon);
%     st2t(19,1:jax2)=imu.*sin(imd*rdcon);
%     st3=zeros(20,jk);
% end;

   
lrgc=cos(lrg*rdcon);
lrgs=sin(lrg*rdcon);


%ship data rate for the cruise is 28 per minute, BUT we interpolated the
%data to have 30 data / min!
for j=1:jk
   jl=30*(j-1)+1;
   jj=min(jax2-1,30*j);
   if jj<jl,
      jl=jj;
   end
%   st3(:,j)=median(st2(:,jl:jj)')';
   st3(:,j)=median(st2t(:,jl:jj)')';
   lrgcm(j)=median(lrgc(jl:jj));
   lrgsm(j)=median(lrgs(jl:jj));
   test(j)=median(sog(jl:jj));
end;

sogx=st3(9,:);
sogy=st3(10,:);

% if ddd>=75
    imux=st3(19,:);
    imuy=st3(20,:);
% elseif ddd<75
%     imux=st3(18,:);
%     imuy=st3(19,:);
% end;

sogm=sqrt(sogx.*sogx+sogy.*sogy)*.514;				%1-min sog, m/s
cogm=mod(atan2(sogy,sogx)/rdcon+360,360);
imum=sqrt(imux.*imux+imuy.*imuy)*.514;				%1-min imet speed, m/s
imdm=mod(atan2(imuy,imux)/rdcon+360,360);
lrgm=mod(atan2(lrgsm,lrgcm)/rdcon+360,360);		%laser ring gyro

% gpstim=st3(2,:);

lat=st3(5,:);
bb=floor(lat/100);
glat=(lat-bb*100)/60;
lat=glat+bb;
clear glat bb

lon=st3(7,:);
bb=floor(lon/100);
glon=(lon-bb*100)/60;
lon=glon+bb;
clear glon bb

tsg=st3(11,:);
tss=st3(12,:);
imrh=st3(13,:);
imta=st3(14,:);
imsol=st3(15,:);
imrain=st3(16,:);
odec=ones(1,length(imrain))*NaN;                 %in m/s

% if ddd>=75
    imir=st3(21,:);
% elseif ddd<75
%     imir=st3(20,:);
% end;


for ii=2:length(lrg),
   if lrg(ii)==0,
      lrg(ii)=lrg(ii-1);
   end;
end

shead=(1:length(tmotion)).*0;
chead=(1:length(tmotion)).*0;

shead(1)=sin(lrg(1)./180.*pi);     %can use also unwrap function to avoid heading problem between 360 and 0 deg.
chead(1)=cos(lrg(1)./180.*pi);

qq=2;

for pp=2:length(tgyroref)-1;
   while tgyroref(pp) >= tmotion(qq),
%      shead(qq)=(sin(lrg(pp+1)./180.*pi)-sin(lrg(pp)./180.*pi))/(tgyroref(pp+1)-tgyroref(pp))*(tgyroref(pp)-tmotion(qq))+sin(lrg(pp)./180.*pi);
%      chead(qq)=(cos(lrg(pp+1)./180.*pi)-cos(lrg(pp)./180.*pi))/(tgyroref(pp+1)-tgyroref(pp))*(tgyroref(pp)-tmotion(qq))+cos(lrg(pp)./180.*pi);
      shead(qq)=(sin(lrg(pp-1)./180.*pi)+sin(lrg(pp)./180.*pi)+sin(lrg(pp+1)./180.*pi))./3;
      chead(qq)=(cos(lrg(pp-1)./180.*pi)+cos(lrg(pp)./180.*pi)+cos(lrg(pp+1)./180.*pi))./3;
      qq=qq+1;
   end;
end;

for ii=2:length(chead),
   if chead(ii)==0,
      chead(ii)=chead(ii-1);
   end;
   if shead(ii)==0,
      shead(ii)=shead(ii-1);
   end;
end;

heading=mod(atan2(shead,chead).*180./pi+360,360);
% ii=find(isnan(heading));heading(ii)=heading(ii(length(ii))+1);

odec2=interp(odec,30);
odec=odec2;clear odec2

sodecx=odec.*cos(lrg./180.*pi);
sodecy=odec.*sin(lrg./180.*pi);
if clearit
clear ii qq pp tmotion tgyro
clear x vul uul tzilch sul st3 st2 st1 sogy sogx sig_sb Tul Tll Rtul Rtll Rsul Rsll
clear Rlul Rlll rdcon prtit prt_it plotit nr1 npz np nl1 lrgsm lrgs lrgcm lrgc jl jk jj 
clear jd jaz jaxx jax jam jade jad j_st j_en j imuy imux hul hr hll g flist imu imd
clear f eps_w dfl dmp2 dmp blk backchk ans
end;