function rainbow(ptlimit)              %last updated  5/17/01
%RAINBOW       A Monte Carlo Simulation of a Rainbow using seawater
%
%    Imagine a large number of of parallel rays hitting a spherical 
%    rain drop. For each ray we plot a single point of light in accordance
%    with the laws of reflection and refraction. For each ray of light we
%    randomly choose a color that will be seen by an observer once it has
%    been refracted & reflected either resulting in a point of light in 
%    the primary or secondary rainbow. Since we are randomly assigning a
%    color from the spectrum to a ray, the process is using a Monte Carlo
%    simulation. In order to get the angles correct we use the refraction
%    index of the colored light for the mediums air and seawater. This also
%    determines the placement of the corresponding dot of light in our picture.   
%
%    The crucial line of code is the expression for N. By changing this
%    line we can simulate rainbows in various mediums.
%
%INPUT:  ptlimit = number of point (rays) to use in the simulation.
%        (If no input is specified 10,000 are used.)
%
%OUTPUT: primary and secondary rainbow arcs
%
%     Original source: "Monte Carlo Computer Simulation of a Rainbow", by
%     Donald Olsen, et al., The Physics Teacher, April 1990, pp. 226-227
%
%     This was the result of an NSF summer program involving eighth and ninth
%     grade students. The original program was written in Applesoft Basic.
%
%
%    David R. Hill, Mathematics Dept., Temple Univ.
%    Philadelphia, PA. 19122    Email: dhill001@temple.edu
%  
%=======================================================================
%CONSTRUCTED for use with NSF Project DEMOs with POSITIVE IMPACT
%=======================================================================

np=0; %counter for number of points plotted

plotit='theta=abs(theta);if theta<= 60,xp=1+1*(theta/60)*(x/b);';
plotit=[plotit 'yp=1*(theta/60)*abs(y/b);'];
plotit=[plotit 'if C==0,pcolor=''r'';end;'];
plotit=[plotit 'if C==1,pcolor=''g'';end;'];
plotit=[plotit 'if C==2,pcolor=''b'';end;'];
plotit=[plotit 'text(xp,yp,''.'',''fontsize'',8,''color'',pcolor,''erasemode'',''none'');'];
%plotit=[plotit 'plot(xp,yp,[''.'' pcolor],''erasemode'',''none'');'];
plotit=[plotit 'drawnow;np=np+1;end;'];

if nargin~=1,ptlimit=10000;end
ptlimit=fix(abs(ptlimit));
%seeding the random number generator
rand('state',sum(100*clock))
%
%the graphics
bkgr='white';
figure('Color',bkgr);
axis([0 2 0 1]),axis(axis);
%axis('square');
axis('off');
hold on;
%loop for generating the rainbow
while np<ptlimit
   b=99;
   while b>=1
      x=-1+2*rand(1);
      y=-1+2*rand(1);
      b=sqrt(x*x+y*y);
   end % while
   %angles of incidence, refraction, & observation
   C=fix(3*rand(1));
   N=1.34 + .01*C;  %IMPORTANT line of code; set to seawater refraction
   I=atan(b/sqrt(1-b*b));
   R=atan(b/sqrt(N*N-b*b));
   T1=(4*R-2*I)*(180/pi);
   T2=(6*R-2*I)*(180/pi)-180;
   %INSTENSITY FACTORS
   RS=(sin(I-R)/sin(I+R))^2;
   RP=(tan(I-R)/tan(I+R))^2;
   I1=(RS*(1-RS)*(1-RS)+RP*(1-RP)*(1-RP))/2;
   I2=(RS*RS*(1-RS)*(1-RS)+RP*RP*(1-RP)*(1-RP))/2;
   if I1>.04*rand(1)
      theta=T1;
      eval(plotit);
   end
   if I2>.02*rand(1)
      theta=T2;
      eval(plotit);
   end
end %while loop for # of points
hold off
disp('RAINBOW is OVER!')