% am 7/25/2012
% reads in omni data
% hardwired filename '/hao/aim3/dmsp/omni/omni_5min1995.asc'
% used for first 90 days therefore can hardwire also number of lines read
% number of lines till '2005  90 23 55' inclusive #of lines 25920
% linear interpolates over fill values
% does the propagation to 1 AU
%
close all

% this part 1 can go outside the loop of the dmsp data, e.g. before the loop
enddoy=90;                 % always reads from doy=1
noflines = enddoy*24*12;   % 12 = 60min/h / 5 min

time = zeros(noflines,5);  % year,doy,ut,min
imf  = zeros(noflines,4);  % bt,bx,by,bz

% specify file to read in & open
file_omni = sprintf('/hao/aim3/dmsp/omni/omni_5min2005.asc');
fid_omni = fopen(file_omni,'r');

% read on line at a time
% checked values with omni website http://omniweb.gsfc.nasa.gov/cgi/nx1.cgi
for nline=1:noflines
   line_in = fgets(fid_omni);
   data = sscanf(line_in,'%d %d %d %d %*f %*f %*f %*f %*f %*f %*f %*f %*f %f %f %*f %*f %f %f');
   %data
   time(nline,1) = data(1);  %  year
   time(nline,2) = data(2);  %  doy
   time(nline,3) = data(3);  %  ut
   time(nline,4) = data(4);  %  min
   time(nline,5) = data(3)*3600+data(4)*60;  %  ut in sec
   imf(nline,1)  = data(5);  %  bt
   imf(nline,2)  = data(6);  %  bx
   imf(nline,3)  = data(7);  %  by
   imf(nline,4)  = data(8);  %  bz
end   % end reading in data from omni file

% find time with fill value
for nline=1:noflines
   if(imf(nline,:)>1000)
     ifound = 0;
     i=0;
     while ifound<1 %find first real value after this fill value 
       if(imf(nline+i,:)<1000)
         ifound = 1;
	 inext=nline+i;
       end
       i=i+1;  % increase by 1
     end
     x1= time(nline-1,2)*86400+time(nline-1,5);  % 24*60*60 sec/day
     x2= time(inext,2)*86400+time(inext,5);  % 24*60*60 sec/day
     x3= time(nline,2)*86400+time(nline,5);
     x= [x1,x2]; 
     for nimf=1:4 
        y= [imf(nline-1,nimf),imf(inext,nimf)];   
        imf(nline,nimf) = interp1(x,y,x3);
     end	
   end
end   % end reading in data from omni file

% close omni file
fclose(fid_omni) 

% end of part1

%part 2 needs to go in the loop over the dmsp data and
% make  idoy_dmsp, utsec_dmsp variable depending on your data

idoy_dmsp = 9              % day number of DMSP data
utsec_dmsp = 45541         % UT in sec of DMSP data

% here the time of the dmsp data is used
% needs to propagate the solar wind measurement from 1 AU to Earth
% calculate index to use for OMNI data based on DMSP data OMNI is 5 min data
% 288: # of index for one day based on 5 min data
% 300: to convert utsec in 5min slices-> 1min/60s/5min=1/300s
i5m = (idoy_dmsp-1)*288 + cast((utsec_dmsp/300.),'int32') + 1;
if (i5m < 5) 
      
%  Fix this later for first 25 min of day
  bx15 = imf(i5m,2);
  by15 = imf(i5m,3);
  bz15 = imf(i5m,4);
  bt15 = imf(i5m,1);
else
  bx15 = (imf(i5m-2,2) + imf(i5m-3,2) + imf(i5m-4,2))/3.;
  by15 = (imf(i5m-2,3) + imf(i5m-3,3) + imf(i5m-4,3))/3.;
  bz15 = (imf(i5m-2,4) + imf(i5m-3,4) + imf(i5m-4,4))/3.;
  bt15 = (imf(i5m-2,1) + imf(i5m-3,1) + imf(i5m-4,1))/3.;
end;
