clear
load data.dat
lambda=[550,625];
day0=data(:,1);
month=data(:,2);
year=data(:,3);
lat=data(:,7);
lon=data(:,6);
day=day2jd(year,month,day0);
hour=data(:,4)+data(:,5)/60;  % GMT time
[nrec,m1]=size(data);
for i=1:nrec
 [az(i),el(i)]=sunflat(year(i),day(i),hour(i),lat(i),lon(i));
 m(i)=airmass(el(i));
end
m=m';
V1=data(:,8);
V2=data(:,9);
ind=find(m<8);
plot(m(ind),log(V1(ind)),'s')
hold on
plot(m(ind),log(V2(ind)),'ro')

A1=polyfit(m(ind),log(V1(ind)),1);
A2=polyfit(m(ind),log(V2(ind)),1);

p=1004;
tau_r=rayleigh(lambda,p);
tau(1)=-A1(1)-tau_r(1);
tau(2)=-A2(1)-tau_r(2);

corr=sundist(day(1));
V01=exp(A1(2))/corr
V02=exp(A2(2))/corr
return
% 
 corr=log(sundist(day(1)));
 tau_t=(-log(V)+log(V0)+corr)./m;
 aot=tau_t-tau_r
