%   'e65/whmatlabtool/datagap/findgap2avehr3.m', also in 'e66/LIshore/*')
%  3/25/2001,  revised from findgap2avehr.m to fill gaps with 'NaN'
% This program 1) find data gaps by converting "calender time" to "sequential time" 
%              2) make hourly average values from 3-min frequency data sampling.
%                "data may consist of missing data point".
%               % convert to hourly data, and fill gaps with 'NaN'.
% note:    this program does 1) fix incorrect data point such as "99.99".
%                            2) fill gap with NaN and cobvert to hourly data 
% by Whuang, 3/27/2001. Test O.K. ( test case using LIshore data )
%

% Input Control ( change the following input control each time ) 

inputfile='p720005.txt'      % input file name; 
y=load(inputfile);            % assign the input data to a variable;
outfile='p720005f.txt'        % give a name to a output file
days=31;                      % total days for the period
month=5;                      % month: 1-January, 2 - Feb, 3- March ...

% -----  You do not need to change following -----

dt=3;                         % sampling time interval inn the input data file;
ymiss= NaN;                   % missing data value; (DO NOT USE 99.99)


% ----  PART 1: FIX INCORRECT DATA
% part 1: fix incorrect data ---
%. assign constrains, and new data set
el=y(:,7);
elmax=8;
eln=el;
nlen=length(el);

% 3. search incorrect data by using the constrain.
ngap=0;
for i=2:nlen;
   
   if ngap==0 & abs(el(i))<elmax;
      eln(i)=el(i);
   end
   
   if abs(el(i))>elmax
      ngap=ngap+1;      
   end
   
   % for only one data point: test O.K;
   if ngap==1 & abs(el(i))<elmax
      eln(i-1)=(el(i-2)+el(i))/2;   
      ngap=0;
   end
   
 %  % for gap >= 2 data points: test O.K 
   if ngap>=2 & abs(el(i))<elmax
       y0=el(i-ngap-1); yend=el(i);
      gradient=(yend-y0)/(ngap+1);
      for j=i-ngap:i-1
         eln(j)= y0+gradient*(j-(i-ngap)+1);
      end       
      ngap=0;   
   end
   
 % end of loop
end
subplot(111)
subplot(311)
plot(el)
xlabel('number of data')
title('e65/whmatlabtool/datagap/fix2all3.m')
legend(inputfile,'3-min raw data') 

'press ENTER to continue!'
pause;

%  ---  PART 2: CONVERT TO HOURLY DATA AND SHOW GAPS 
ntothr=days*24;
mmddyyyy=y(:,2);
hhmmss=y(:,3);
%el=y(:,7);

% 2. constrains: data interval =3 min;
% time label only track month-day-hour to convert hourly-average value
mmdd=mmddyyyy/10000;
hh=hhmmss/10000;

% 3) open temp output file to check output.
fid = fopen('temp.txt','w'); % 
len=length(y);
for i= 1:len
fprintf(fid,'%10.0f %10.0f %10.4f \n', mmdd(i), hh(i),eln(i));
end
fclose(fid)

%pause
% 4) reload from temp file; find hourly average value, and write to a file:
%clear;
load temp.txt;
mmdd=temp(:,1);
hh=temp(:,2);
mmddhh=(mmdd-month*100)+hh/24;

el=temp(:,3);
len=length(el);

%fid = fopen('testhrly.txt','w'); % 
fid = fopen('temp1.txt','w'); % 
nhr=0;
vhr=mmddhh(1);
t=mmddhh*0;
y=t*0;

ncount=0;
for i= 2:len-1
   if mmddhh(i)==vhr;
      nhr=nhr+1;
   end
   
   if mmddhh(i+1)>vhr;
      elave=mean(el(i-nhr+1:i));
      ddhhave=mean(mmddhh(i-nhr+1:i))-1;
      fprintf(fid,'%10.0f %10.0f %10.2f %10.2f \n', mmdd(i), hh(i),ddhhave, elave);
      vhr=mmddhh(i+1); 
      nhr=0;
      ncount=ncount+1;
       end
end

 fclose(fid)
% 4'. re-load, and rewrite to include NaN value for missing data (large gaps)

yt=load('temp1.txt');
fid1 = fopen(outfile,'w'); 

z=zeros([31*24,4]);
z=z+ymiss;
nct=1;

for i=1:31
   ii=i+500;
   for j=1:24;
      jj=j-1;
      for k=1:length(yt);
         if yt(k,1)==ii & yt(k,2)==jj;
            z(nct,1)=ii; z(nct,2)=jj; z(nct,3)=yt(k,3); z(nct,4)=yt(k,4);
         end
      end
       nct=nct+1;
   end
 end

   
 for i=1:ntothr;
    fprintf(fid1,'%10.0f %10.0f %10.4f %10.2f \n', z(i,1), z(i,2),z(i,3), z(i,4));  
 end
 
 
 fclose(fid1);

% 

 
 % 5. plot to compare 
  
 subplot(312)
 plot(mmddhh,el)
 xlabel('Time (days from 5/1/2000)')
  legend('fix small gap, show large gaps in raw data')

 
 subplot(313)
  s=load(outfile);
  plot(s(:,3),s(:,4))
 xlabel('Time (days from 5/1/2000)')
 legend( outfile,'hourly data, show gaps')

 






