function sun(zen)
%clear
%zen=87.5;
dsun=0.5;

[x,y]=meshgrid(zen-dsun/2:0.002:zen+dsun/2,zen-dsun/2:0.002:zen+dsun/2);

el=90-y;
del=refraction(el);
el=el+del;   
y=90-el;

lambda=[600,550,350];
AOT=0.03;
alfa=1;
AOT=AOT*(lambda/500).^-alfa;
m=airmass(90-y);
tau_r=rayleigh(lambda,1000);
load solar.dat
L=10000./solar(:,1);
I0new=10*solar(:,2)*10000./L.^2;
L=L*1000;
I0=interp1(L,I0new,lambda);
%I0=[0.6,1,0.6];
for i=1:length(lambda)
    I(:,:,i)=I0(i).*sin(el*pi/180).*exp(-m*(tau_r(i)+AOT(i)));
end 

MAX=max(max(max(I)));
I=I/MAX;
zen_cor=zen-refraction(90-zen);
zen_min=min(min(zen));
zen_max=max(max(zen));
dzen_cor=max(max(y))-min(min(y));
for i=1:length(x)
    for j=1:length(y);
        %dist=sqrt((x(i,j)-zen)^2+(y(i,j)-zen_cor)^2);
        dist=((x(i,j)-zen)/(dsun/2))^2+((y(i,j)-zen_cor)/(dzen_cor/2))^2;
        if dist>1    %dsun/2
            I(i,j,:)=[0,0,1];
        end
    end
end  
%h=polar([0 2*pi], [0 85],'b');
%delete(h);
%hold on;
%subplot(2,2,3)
imshow(I)
h=title([num2str(zen)])

set(gcf,'Color',[0,0,1])
return
Fi=[0:1:360]*pi/180;
r1=dsun*ones(size(Fi));
[X,Y]=pol2cart(Fi,r1);
fill(X+zen,Y+zen,'b')