function rainbow(ptlimit)              %last updated  6/13/01
%RAINBOWFULL       A Monte Carlo Simulation of a Rainbow 
%
%    In  this version 7 colors are used.
%
%    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 water. 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 the
%    refraction index values 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
%=======================================================================

%REFRACTION INDEX DATA
refracindex=[1.33 1.34 1.35 1.36 1.37 1.38 1.39];
%the preceding data are "fake"; they is used to be able to distinguish the color bands
%===================================================================================
%refracindex=[1.3309 1.3328 1.3339 1.3349 1.3392 1.3419 1.3431];
%the preceding data are close to the actual refraction indexes
%using these values it is difficult to distinguish the color bands
%===================================================================================
%NOTE: index in ROY G BIV order

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=[1 0 0];end;'];%RED
plotit=[plotit 'if C==1,pcolor=[1 .5 0];end;'];%ORGANGE
plotit=[plotit 'if C==2,pcolor=[1 1 0];end;'];%YELLOW
plotit=[plotit 'if C==3,pcolor=[0 1 0];end;'];%GREEN
plotit=[plotit 'if C==4,pcolor=[0 0 1];end;'];%BLUE
plotit=[plotit 'if C==5,pcolor=[.5 0 .5];end;'];%INDIGO
plotit=[plotit 'if C==6,pcolor=[.67 0 1];end;'];%VOILET
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 %code to make sure for pair (x,y) 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(7*rand(1));%modified to use seven colors
   %N=1.33 + .01*C;  %IMPORTANT code; sets refraction idexes; 
                     % original statement for red, green, blue rainbow
   N=refracindex(C+1);%IMPORTANT code; sets refraction idexes; +1 because used as subscript
   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('RAINBOWFULL is OVER!')