%------ read_shp_HiWinGS_2013.m

%	reads daily SCS ship data system files for HiWinGS / RV Knorr

%	These variables must exist in the current workspace
%		hhh: numerical value of current hour (00-23)
%		ddd: numerical value of current day of year (0-365)
%		yyyy: string variable of current year (e.g.'2014')
%		way_raw_data_flux: string path to data folder
%		jd_pc: 1 min timestamp from slow met data (read_hp_HiWinGS_2013.m)
%       clearit: sets clear flag

%	Script creates these arrays:
%		gpstime, lat, lon, sogx, sogy, tsg, tss, imrh, imta, imsol, imrain
%		odecm, imux, imuy, imir, sogm, cogm, imum, imdm, lrgm, jd_scs, lrg

% File Format:  approx 0.1Hz
% IDX     TAG    Date        Time_GMT  Dec_Lat    Dec_Lon     SPD     HDT    COG    SOG   WXTP_Ta WXTS_Ta WXTP_Pa WXTS_Pa WXTP_Ri WXTS_Ri WXTP_Rc WXTS_Rc WXTP_Dm WXTS_Dm WXTP_Sm WXTS_Sm WXTP_Ua WXTS_Ua WXTP_TS WXTS_TS WXTP_TD WXTS_TD SWR    LWR   ND   SBE45S  SBE48T  FLR       SSV      Depth12 Depth35 
% 3031875 $dsLog 2013/10/15  00:30:31  58.380150  -45.270833  000.09  052.7  333.2  1.53  7.6     7.4     1009.3  1009.2  0.0     0.0     13.07   16.27   004     359     19.6    19.4    81.7    84.4    19.53   19.26   58.99   53.99   15.652 345.6 NaN  34.9598 7.0424  00095.50  1480.41  00.0    00.00
% 1 2 3          4    5  6   7  8  9   10         11          12      13     14     15   16       17      18      19      20      21      22      23      24      25      26      27      28      29      30      31      32      33      34     35    36   37      38      39         40      41      42

%% initialize
disp('read_shp_hiwings2013.m')
fclose('all');
jd=sprintf('%03i',ddd); % day of year string
hr=sprintf('%02i',hhh); % hour string

nr2=43; % number of variables

% Full path & file name string for Windows computer
% dfl=[way_raw_data_flux, yyyy(3:4),jd,'\mot0' yyyy(3:4) jd hr '_raw.txt'];

% Full path & file name string for Mac
dfl=[way_raw_data_flux, yyyy(3:4),jd,'/scs0', yyyy(3:4) jd hr '_raw.txt'];

disp(['Reading scs file from doy ',jd,'  hour ',hr]);
flist=fopen(dfl,'r');

temp={};    % empty cell array for textscan
stx=[];     % empty array for output

%% read raw text files
if flist>0,             % if the file exists,
    while ~feof(flist)  % read file into cell array
        try
        	% temp is 1x(nr2) cell array of ~360x1 arrays
            temp=textscan(flist,'%2f%2f%3f %*6c %4f/%2f/%2f %2f:%2f:%2f %f %f %f %f %f %f %f %f %f %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)'];	% 42x~360 array
        catch
            for ii=1:nr2    % length of last cell reflects missing values from any/all fields
                 if length(temp{1,ii})~=length(temp{1,nr2})
                     temp{1,ii}(lnr2+1)=[]; %truncate length of all cells to last cell
                 end;
            end;
            if isempty(temp{1,1}); continue; end;   % continue if no data yet
            stx=[stx cell2mat(temp)'];
        end;
    end;
    fclose(flist);
else  % no data default
    stx=ones(nr2,360)*NaN;
end;

%% variable definitions
% jd_scs=datenum(stx(6,:),stx(5,:),stx(4,:),stx(7,:),stx(8,:),stx(9,:))-datenum(Vdate(1),1,0);   %Day-Of-Year
% Dec_Lat=stx(10,:);  %Latitude	deg	C-Nav GPS receiver
% Dec_Lon=stx(11,:);  %Longitude	deg	C-Nav GPS receiver
% SPD=stx(12,:); 
% HDT=stx(13,:);      %heading
% COG=stx(14,:);      %Course Over Ground	deg	C-Nav GPS receiver
% SOG=stx(15,:);      %Speed Over Ground	knots	C-Nav GPS receiver
% WXTP_Ta=stx(16,:);  %Air Temp	C	Port VaisalaWXT520
% WXTS_Ta=stx(17,:);  %Air Temp	C	Starboard VaisalaWXT520
% WXTP_Pa=stx(18,:);  %Pressure	hPa	Port VaisalaWXT520, uncorrected for height
% WXTS_Pa=stx(19,:);  %Pressure	hPa	Starboard VaisalaWXT520, uncorrected for height
% WXTP_Ri=stx(20,:);  %Rainfall Intensity	mm/h	Port VaisalaWXT520
% WXTS_Ri=stx(21,:);  %Rainfall Intensity	mm/h	Starboard VaisalaWXT520
% WXTP_Rc=stx(22,:);  %Cumulative Rainfall	mm	Port VaisalaWXT520
% WXTS_Rc=stx(23,:);  %Cumulative Rainfall	mm	Starboard VaisalaWXT520
% WXTP_Dm=stx(24,:);  %Relative Wind Direction	deg	Port VaisalaWXT520
% WXTS_Dm=stx(25,:);  %Relative Wind Direction	deg	Starboard VaisalaWXT520
% WXTP_Sm=stx(26,:);  %Relative Wind Speed	m/s	Port VaisalaWXT520
% WXTS_Sm=stx(27,:);  %Relative Wind Speed	m/s	Starboard VaisalaWXT520
% WXTP_Ua=stx(28,:);  %Relative Humidity	%	Port VaisalaWXT520
% WXTS_Ua=stx(29,:);  %Relative Humidity	%	Starboard VaisalaWXT520
% WXTP_TS=stx(30,:);  %True Wind Speed	m/s	Port VaisalaWXT520, corrected for ship's gyro, SOG, and COG
% WXTS_TS=stx(31,:);  %True Wind Speed	m/s	Starboard VaisalaWXT520, corrected for ship's gyro, SOG, and COG
% WXTP_TD=stx(32,:);  %True Wind Direction	deg	port VaisalaWXT520, corrected for ship's gyro, SOG, and COG
% WXTS_TD=stx(33,:);  %True Wind Direction	deg	Starboard VaisalaWXT520, corrected for ship's gyro, SOG, and COG
% SWR=stx(34,:);      %Short Wave Radiation	W/m2	IMET (Eppley Laboratories, Inc. pyronometer)
% LWR=stx(35,:);      %Long Wave Radiation	W/m2	IMET (Eppley Laboratories, Inc. pyronometer)
% SBE45S=stx(36,:);   %Salinity ppth SBE45
% SBE48T=stx(37,:);   %Sea Surface Temperature	C	SeaBird SBE48
% FLR=stx(39,:);      %Fluorescence	mV	WetLabs Wet-StarWS3S-901P
% SSV=stx(40,:);      %Surface Sound Velocity	m/s	Calculated from SSTMP and SSCND
% Depth12=stx(41,:);
% Depth35=stx(42,:);

%% select desired output variables
% stx(:,any(isnan(stx),1))=[]; % remove lines containing NaNs (bad lines, too long...)

lat =     stx(10,:);        % Latitude	deg	C-Nav GPS receiver
lon =     stx(11,:);        % Longitude	deg	C-Nav GPS receiver
cog =     stx(14,:);        % Course Over Ground	deg	C-Nav GPS receiver
sog =     stx(15,:)*0.514;  % Speed Over Ground	knots	C-Nav GPS receiver   %converted from knot to m/s
tsg =     stx(36,:);        % Salinity ppth SBE45
tss =     stx(37,:);        % Sea Surface Temperature	C	SeaBird SBE48
rh =      stx(29,:);        % Relative Humidity	%	Starboard VaisalaWXT520
ta =      stx(17,:);        % Air Temp	C	Starboard VaisalaWXT520
solar =   stx(34,:);        % Short Wave Radiation	W/m2	IMET (Eppley Laboratories, Inc. pyronometer)
imrain =  stx(21,:);        % Rainfall Intensity	mm/h	Starboard VaisalaWXT520
lrg =     stx(13,:);        % heading
odec =    stx(15,:)*NaN;    % speedlog - not working on HiWinGS
imu =     stx(30,:);        % True Wind Speed	m/s	Port VaisalaWXT520, corrected for ship's gyro, SOG, and COG
imd =     stx(32,:);        % True Wind Direction	deg	port VaisalaWXT520, corrected for ship's gyro, SOG, and COG
ir  =     stx(35,:);        % Long Wave Radiation	W/m2	IMET (Eppley Laboratories, Inc. pyronometer)
impres =  stx(18,:);        % Pressure	hPa	Port VaisalaWXT520, uncorrected for height
imur   =  stx(26,:);        % Relative Wind Speed	m/s	Port VaisalaWXT520
imdr   =  stx(24,:);        % Relative Wind Direction	deg	Port VaisalaWXT520

gpstim=hhh+stx(1,:)./60+(stx(2,:)+stx(3,:)/1000)./60./60;      %Decimal time in hours

for iii=2:length(gpstim),	% Fix bad timestamps
    if gpstim(iii)<=gpstim(1),
        gpstim(iii)=gpstim(iii-1)+gpstim(2)-gpstim(1);
    end;
end;

%% find zero values of rh, sal, sst, lat, lon and ir - set to prior value
ii=find(rh==0);rh(ii)=rh(ii-1);
ii=find(tsg==0);tsg(ii)=tsg(ii-1);
ii=find(tss==0);tss(ii)=tss(ii-1);
ii=find(lat==0);lat(ii)=lat(ii-1);
ii=find(lon==0);lon(ii)=lon(ii-1);
ii=find(ir==0);ir(ii)=ir(ii-1);
% find bad data flags in wind and set to NaN
ii=find(imu==-99);imu(ii)=NaN;ii=find(imur==-99);imur(ii)=NaN;
ii=find(imd==-99);imd(ii)=NaN;ii=find(imdr==-99);imdr(ii)=NaN;

%% Format selected variables into new array of ~10 sec data
st2=zeros(16,length(gpstim));	 % 16x~360 array
st2(1,:)=gpstim;
st2(2,:)=lat;
st2(3,:)=lon;
st2(4,:)=cog;
st2(5,:)=sog;
st2(6,:)=tsg;
st2(7,:)=tss;
st2(8,:)=rh;
st2(9,:)=ta;
st2(10,:)=solar;
st2(11,:)=imrain;
st2(12,:)=lrg;
st2(13,:)=odec;
st2(14,:)=imu;
st2(15,:)=imd;
st2(16,:)=ir;

st2=st2';	% rotate to ~360x16 array

%% convert ship speed, wind speed and heading in st2 to N E component vectors
rdcon=pi/180;
st2(:,4)=sog'.*cos(cog'*rdcon);    	% SOG N vector component
st2(:,5)=sog'.*sin(cog'*rdcon);    	% SOG E vector component
st2(:,14)=imu'.*cos(imd'*rdcon);   	% TWSPD N vector component
st2(:,15)=imu'.*sin(imd'*rdcon);   	% TWSPD E vector component
lrgc=cos(lrg'*rdcon);        		% HDG N vector component
lrgs=sin(lrg'*rdcon);        		% HDG E vector component

if clearit
    clear nr2 temp jd hr dfl
    clear stx ii iii lat lon sog cog tsg tss rh ta solar imrain odec
    clear imu imd ir impres imur imdr
end;

%% average variables to 1 min timebase
% jd_ref=jd_pc;						% use slow met time stamp (dec DOY @ 1 min)
jd_ref = ddd + hhh/24 + (0:1:59)/1440; % use exact 1-min timestamp
jd_scs = gpstim./24+ddd;				% convert ship timestamp to dec DOY
st3 = zeros(length(jd_ref),16);		% 60x16 empty array for 1 min data
lrgsm = zeros(1,length(jd_ref));      % empty arrays for 1 min heading
lrgcm = zeros(1,length(jd_ref));

% Average or interp data to jd_ref
for ii=1:16
%	st3(:,ii)=interp1(jd_scs(~isnan(jd_scs)),st2(~isnan(jd_scs),ii),jd_ref);
	st3(:,ii)=intavg(jd_scs(~isnan(jd_scs)),st2(~isnan(jd_scs),ii),jd_ref);
end
% lrgcm=intavg(jd_scs,lrgc,jd_ref); 
% lrgsm=intavg(jd_scs,lrgs,jd_ref); 
lrgcm=interp1(jd_scs(~isnan(jd_scs)),lrgc(~isnan(jd_scs)),jd_ref); 
lrgsm=interp1(jd_scs(~isnan(jd_scs)),lrgs(~isnan(jd_scs)),jd_ref); 

if length(jd_ref)<60	% pad to 60 min with NsN if necessary
    st3(length(jd_ref):60,:)=NaN;
    lrgcm(length(jd_ref):60)=NaN;
    lrgsm(length(jd_ref):60)=NaN;
end

% varibles synchronized with PSD slow met data
% all arrays are 1x60
gpstime=st3(:,1)';      % decimal hour corresponding to jd_ref
lat=st3(:,2)';%./100;   % lat
lon=st3(:,3)';%./100;   % lon
sogx=st3(:,4)';         % SOG N
sogy=st3(:,5)';         % SOG E
tsg=st3(:,6)';          % Salinity
tss=st3(:,7)';          % Tsea
imrh=st3(:,8)';         % RH
imta=st3(:,9)';         % Tair
imsol=st3(:,10)';       % solar
imrain=st3(:,11)';      % Rain
odecm=st3(:,13)';       % Speedlog
imux=st3(:,14)';        % Twspd N
imuy=st3(:,15)';        % Twspd E
imir=st3(:,16)';        % Longwave
sogm=sqrt(sogx.*sogx+sogy.*sogy);               % computed 1-min sog
cogm=mod(atan2(sogy,sogx)/rdcon+360,360);       % computed 1-min sog
imum=sqrt(imux.*imux+imuy.*imuy);               % computed 1-min Twspd, m/s
imdm=mod(atan2(imuy,imux)/rdcon+360,360);       % computed 1-min Twdir, deg
lrgm=mod(atan2(lrgsm,lrgcm)/rdcon+360,360);     % computed heading (laser ring gyro)

if clearit
    clear st2 st3 lrgc lrgs gpstim jd_ref
end;
