% Read DMS files for HiWinGS 2013 data

function [jd_UH,AccX_UH,AccY_UH,AccZ_UH,RateX_UH,RateY_UH,RateZ_UH,U_UH,V_UH,W_UH,Tsonic_UH,LagPulse,Xlag]=read_dms_HiWinGS_2013(ddd,hhh,yyyy,varPSD)

% disp(['read_DMS: ', cruise]);
 jd=sprintf('%03i',ddd);
hr=sprintf('%02i',hhh);
fclose('all');

%declare arrays
jax=1;                              %line counter
nr2=13;
dmsdata=ones(nr2+1,36000)*NaN;


% for jam=0:23,					    %cycle thru 24 hourly stats files
    
    e=['G:\Data\CRUISES\HiWinGS_2013\Knorr\UH_DMS\' jd '\DMS_', yyyy '_' jd '_' hr '00.txt'];
    %e=[way_raw_data_flux, 'day',jd,'\met0' year(3:4) jd hr '_raw.txt'];
    disp(['Reading DMS file from hour ',e]);
    flist=fopen(e,'r');
    temp=[];stx=[];
    
    if flist>0,					    %if the file exists,
        while ~feof(flist)
            try
                temp =textscan(flist,'%f %*f %*f %*f %*f %*f %f %f %f %f %f %f %f %f %f %f %f %*f %f %*f %*f %*f %*f %*f %*[^\n]', 'delimiter', ',', 'headerlines', 1,  'emptyvalue', NaN);
                stx=[stx cell2mat(temp)'];
            catch
                [lnr2,cnr2]=size(temp{1,nr2});
                for ll=1:nr2
                     if length(temp{1,ll})~=lnr2
                         temp{1,ll}(lnr2+1)=[];
                     end;
                end;
                if isempty(temp{1,1}); continue; end
                if isempty(temp{1,1});
                    disp('Should never do this!')
                    temp{1,nr2}=temp{1,nr2}';
                end;
             stx=[stx cell2mat(temp)'];
            end;
        end;  %end while
    
        [nr1,nl1]=size(stx);
        hhrs=ones(1,nl1)*hhh;
        dmsdata(:,jax:jax+nl1-1)=[stx' hhrs']';
        jax=jax+nl1;
    
    else
        disp(sprintf('\b<--File not found.'))
        stx=ones(nr2,36000)*NaN;
%         stx(1,:)=ones(1,3600)*jam; stx(2,:)=ones(1,3600)*jam;
        [nr1,nl1]=size(stx);
        hhrs=ones(1,nl1)*jam;
        dmsdata(:,jax:jax+nl1-1)=[stx' hhrs']';
        jax=jax+nl1;
    end;      %     end if list
  

% if ddd==282
%     dmsdata(1,1)= 24278400.0239;
% end
% jd_UH= ddd+(dmsdata(1,:)-dmsdata(1,1))/3600/24; % time in decimal yearday
jd_UH= 1+dmsdata(1,:)/86400; % time in decimal yearday
U_UH=(dmsdata(2,:)-2.5).*12;     % scale wind to physical units
V_UH=(dmsdata(3,:)-2.5).*12;  
W_UH=(dmsdata(4,:)-2.5).*12; 

uv = [U_UH;V_UH];
theta = (-30*(pi/180));
Tm = [(cos(theta)) (-sin(theta));(sin(theta)) (cos(theta))];
UVrot = Tm*uv;
U_UH = UVrot(1,:);
V_UH = UVrot(2,:);
% Reverse sign of U to agree with motion coordinates (forward positive)
U_UH = -U_UH;
   
Spd_Snd=dmsdata(5,:); 
%compute sonic temperature, NOAA R2 formula
SpdSnd = Spd_Snd.*16 + 290;  % speed of sound in m/s
Tsonic_UH = (SpdSnd.*SpdSnd + (U_UH.*U_UH + W_UH.*W_UH)/2. + V_UH.*V_UH)/403-273.15;    

% convert motion to physical units and redefine axes according to
% standard convention.
%accx  : forward is positive 
%accy  : to port is positive
%accz  : up is positive
%ratex : port up is positive phi
%ratey : bow down is positive theta
%ratez : bow to port is positive psi (right-handed)

RateX_UH = dmsdata(7,:)'.*(1000/50.134).*(pi/180);
RateY_UH = dmsdata(6,:)'.*(1000/50.034).*(pi/180);
RateZ_UH = -dmsdata(8,:)'.*(1000/49.975).*(pi/180);
AccX_UH = dmsdata(10,:)'.*(9.80665/7.510);
AccY_UH = dmsdata(9,:)'.*(9.80665/7.492);
AccZ_UH = -dmsdata(11,:)'.*(9.80665/3.747);

% 
% figure;
% subplot(3,2,1);plot(ddd+tref/3600/24,accx,jd_UH,AccX_UH,'r');title('accx');
% subplot(3,2,2);plot(ddd+tref/3600/24,accy,jd_UH,AccY_UH,'r');title('accy');
% subplot(3,2,3);plot(ddd+tref/3600/24,accz,jd_UH,AccZ_UH,'r');title('accz'); 
% subplot(3,2,4);plot(ddd+tref/3600/24,ratex,jd_UH,RateX_UH,'r');title('ratex');
% subplot(3,2,5);plot(ddd+tref/3600/24,ratey,jd_UH,RateY_UH,'r');title('ratey');
% subplot(3,2,6);plot(ddd+tref/3600/24,ratez,jd_UH,RateZ_UH,'r');title('ratez');
if ~isnan(varPSD)
%     Xlag=findlag(varPSD,AccZ_UH,20,0);
    Xlag=findlag(varPSD,W_UH,20,0);
else
    Xlag=0;
end
LagPulse=dmsdata(12,:);
% % BlankPulse=dmsdata(13,:); %related to DMS
% figure;plot(jd_UH,LagPulse,'+');grid;