clear
close all
name=['MICROTOPS2016.dat'];
fid=fopen(name,'r');           % open the data's file
if fid<0
 disp('No file')
 return
end 
str=fgets(fid);                % omission the first line in file
str=fgets(fid);                % get the second line in file

for i=1:5
 x1=str((7+16*(i-1)):(7+16*(i-1)));
 x2=str((9+16*(i-1)):(11+16*(i-1)));
 x3=str((13+16*(i-1)):(15+16*(i-1)));
 x1=str2num(x1);
 x2=str2num(x2);
 x3=str2num(x3);
 LNV0(i)=(x1+x2/1000)*10.^(x3);     % signal calibration constant
end
                                    % get the next line in file
str=fgets(fid);

for i=1:5
 x1=str((4+13*(i-1)):(4+13*(i-1)));
 x2=str((6+13*(i-1)):(8+13*(i-1)));
 x3=str((10+13*(i-1)):(12+13*(i-1)));
 x1=str2num(x1);
 x2=str2num(x2);
 x3=str2num(x3);
 cal(i)=(x1+x2/1000)*10.^(x3);      % calibration constant
end

str=fgets(fid);                     % get the next line in file

for i=1:2
 x1=str((3+12*(i-1)):(3+12*(i-1)));
 x2=str((5+12*(i-1)):(7+12*(i-1)));
 x3=str((9+12*(i-1)):(11+12*(i-1)));
 x1=str2num(x1);
 x2=str2num(x2);
 x3=str2num(x3);
 k0(i)=(x1+x2/1000)*10.^(x3);
end
 K=k0(1);                           % constant K
 B=k0(2);                           % cinstant B
 x1=str(27:27);
 x2=str(29:30);
 x1=str2num(x1);
 x2=str2num(x2);
 C=x1+x2/100;                       % constant C

str=fgets(fid);                     % get the next line in file
 x1=str(7:8);
 x2=str(10:12);
 x3=str(14:16);
 x1=str2num(x1);
 x2=str2num(x2);
 x3=str2num(x3);
 POFFS=(x1-x2/1000)*10.^(x3);       % first pressure constant


 x1=str(25:25);
 x2=str(27:29);
 x3=str(31:33);
 x1=str2num(x1);
 x2=str2num(x2);
 x3=str2num(x3);
 PSCALE=(x1+x2/1000)*10.^(x3);      % second pressure constant

 str=fgets(fid);                     % omission 3 line in file
 str=fgets(fid);
 str=fgets(fid);
 a=strfind(str,'SIG');
 str=str(a(1):a(end)+10);
 a=strfind(str,'SIG');
 b=strfind(str,',');
 for i=1:length(a)
  lambda(i)=str2num(str(a(i)+3:b(i)-1));
 end
 str=fgets(fid);
 i=0;
 while (str~=-1)
    i=i+1;
    x=str;
    serial(i)=str2num(x(1:5));                     % serial number
    year(i)=str2num(x(13:16));                     % year
    month(i)=str2num(x(7:8));                      % month
    day(i)=str2num(x(10:11));                        % day
    day0(i)=day(i);
    hour(i)=str2num(x(18:19));                       % hour
    hour0(i)=hour(i);
    minute(i)=str2num(x(21:22));                     % minute
    sec(i)=str2num(x(24:25));                        % secund
    hour(i)=hour0(i)+1/60*minute(i)+1/3600*sec(i);  % calculate hour in decimal
    aa=findstr(str,',');
%    strend=aa(length(aa));
    str_new=str(27:end);
%    aNaN=findstr(str_new,'Na');
    str_new=strrep(str_new,'Na','NaN');
    str_new=strrep(str_new,'NaNN','NaN');
    
%    if isempty(aNaN)==1
      X=str2num(str_new);
%      keyboard
      lat(i)=X(1);                       % latitude
      lon(i)=X(2);                       % longtitude
      alt(i)=X(3);                       % altitude
      p(i)=X(4);                         % pressure
      zen(i)=X(5);                        % sun azimuth angle
      m(i)=X(6);                         % air mass
      dist(i)=X(7);                      % Earth-sun distance ratio
      t(i)=X(8); 
%      disp(i)                        % temperature of microtops
      sig(i,:)=X(10:14);                    % signal
      std(i,:)=X(15:19);                    % standar deviation of signal
      R380_440(i)=X(20);
      R440_500(i)=X(21);                  % signal 440 to 500 ratio
      R500_670(i)=X(22);                  % signal 500 to 675 ratio
      R675_870(i)=X(23);                  % signal 675 to 870 ratio
      AOT(i,:)=X(24:28);                    % Aerosol optical depth
      PW(i)=X(29);
      if length(X)==30
       flag(i)=X(end);
       else
       flag(i)=NaN;
      end
   str=fgets(fid);
 end
fclose(fid);
 %cd(path_old)
  nrec=length(lat);                     % number of records   
 for i=1:nrec
    julian_day(i)=day2jd(year(i),month(i),day0(i));
 end                                % calculate julian day

doy=julian_day+hour/24;
data.serial=serial(1);
data.year=year;
data.jd=julian_day;
data.month=month;
data.day=day;
data.hour=hour0;
data.minute=minute;
data.sec=sec;
data.time=hour;
data.doy=doy;
data.sig=sig;
data.std=std;
data.AOT=AOT;
data.AOTraw=AOT;
data.lnv0fact=LNV0;
data.calPW=[K,B,C];
data.lambda=lambda;
data.lat=lat;
data.lon=lon;
data.alt=alt;
data.p=p;
data.flag=flag;

for i=1:nrec             % calculate sun position -azimyth and elevation angle
 [az(i),el(i)]=sunpos(data.year(i),data.jd(i),data.time(i),data.lat(i),data.lon(i));
end
m=airmass(el);              % calculate airmass
data.az=az;
data.el=el;
data.m=m;

lnv0=[7.6641    6.8222    7.0034    6.5164    7.0763]; %WRND 2010
if year(i)==205
 lnv0=[7.7505    6.7395    6.7182    7.0762    6.6152];
end 
for i=1:length(lambda)                    % total optical depth
 tau_total(:,i)=(LNV0(i)+log(sundist(data.jd'))-log(eps+data.sig(:,i)))./m';
end
tau_R=rayleigh(data.lambda,data.p');     % Rayleigh optical depth
tau_O3=opt_ozon(lambda,300);

for i=1:nrec
 tau_AOT(i,:)=tau_total(i,:)-tau_R(i,:)-tau_O3;      % calculate Aerosol Optical Depth
end
data.AOT=tau_AOT;
ind1=find(lambda==500);
ind2=find(lambda==870);
if (isempty(ind1)==0)&(isempty(ind2)==0)
 alfa=log(data.AOT(:,ind1)./data.AOT(:,ind2))/log(data.lambda(ind2)/data.lambda(ind1));
else
 alfa=NaN*ones(size(data.time));
end 
for i=1:nrec
   A=polyfit(log(data.lambda),log(data.AOT(i,:)),1); 
   ALFA(i)=-A(1);
end
data.angs500_870=alfa';
data.angs=ALFA;

% Water vapour channel
ind936=find(lambda==936);
if isempty(ind936)==0
 data.AOT(:,ind936)=NaN;
 ind1020=find(lambda==1020);
 ind870=find(lambda==870);
 PWangs=log(data.AOT(:,ind870)./data.AOT(:,ind1020))/log(1020/870);
 Cangs=(1020/936).^PWangs;   % or 1.16 from manual
 
 term1=AOT(:,ind1020).*data.m'.*(1-Cangs)-log(data.sig(:,ind936)./data.sig(:,ind1020)*LNV0(ind1020)/LNV0(ind936));
 term2=data.calPW(1)*data.m.^data.calPW(2);
 data.PWraw=PW;
 data.PW=(term1./term2').^(1/data.calPW(2));
end 

% wybiera dane z serii pomiarowej

time=data.doy;
t0=time(1);
k=0;
dt=0.05/24;
chan=2;
AOT=data.AOT(:,chan);
while 1
   k=k+1; 
   ind=find((time>=t0)&(time<=t0+dt)); 
   if length(ind)>=3
%   w=(1./data.std(ind,chan)).^2;
%   AOT(k)=sum(data.AOTnew(ind,chan).*w)/sum(w); 
   [val(k),index]=min(AOT(ind));
   IND(k)=ind(index);
   else
    k=k-1;
   end
   try
    t0=time(ind(end)+1);
   catch
       break
   end
end    

data.year=data.year(IND);
data.jd=data.jd(IND);
data.month=data.month(IND);
data.day=data.day(IND);
data.hour=data.hour(IND);
data.minute=data.minute(IND);
data.sec=data.sec(IND);
data.time=data.time(IND);
data.doy=data.doy(IND);
data.sig=data.sig(IND,:);
data.std=data.std(IND,:);
data.AOT=data.AOT(IND,:);
data.AOTraw=data.AOTraw(IND,:);
data.lnv0fact=LNV0;
data.lambda=lambda;
data.lat=data.lat(IND);
data.lon=data.lon(IND);
data.alt=data.alt(IND);
data.p=data.p(IND);
data.az=data.az(IND);
data.el=data.el(IND);
data.m=data.m(IND);
data.flag=data.flag(IND);
data.angs500_870=data.angs500_870(IND);
data.angs=data.angs(IND);
if isfield(data,'PW');
 data.PW=data.PW(IND);
 data.PWraw=data.PWraw(IND);
end
save MICROTOPS data
