program photons8;  {Program to simulate laser tweezers by Steven B. Smith.  Compiled using Borland Turbo Pascal}
uses crt, graph;
const
pLength = 10;       {number of pixels depicting photon}
rayNum = 200;       {number of photons on screen}
relInd= 1.58/1.334; {relative index change from polyStyrene/water}
LightUmph = 5;      {nudge factor for light force}

type
path = array[0..pLength,1..2] of real;
photonType = record
             color:byte; {color on screen}
             Px,Py:real; {x and y momentum}
             track:path; {positions of pixels in photon}
             endPtr:byte;{index of leading pixel}
             inHiIndex:boolean; {photon inside high-index medium}
             end;
var
grDriver,grMode,xMax,xMid,yMax,yMid,loopCtr,yOffset,xPos,yPos,Length,
oldWidth,FxPos,FyPos,beadR,photoCount,transmitRate:integer;
dualBeam, colimated,menuOn,refracted,reflected:boolean;
photon:array[-rayNum..rayNum] of photonType;
beadX,beadY,spread,xforce,yforce,angle,targetX,targetY,rate:real;
ch:char;
var reflCtr,refrCtr:integer;


procedure lineBlob (xPos, yPos,length:integer; angle:real);
var xEnd,yEnd:integer;  {makes an "arrow" to show vector}
begin
if length >0 then
  begin
  xEnd:= xPos+round(length*cos(angle));
  yEnd:= yPos+round(length*sin(angle));
  line(xPos,yPos,xEnd,yEnd);
  line(xPos+1,yPos+1,xEnd+1,yEnd+1);
  circle(xEnd,yEnd,3);
  circle(xEnd,yEnd,2);
  end;
end;

procedure showForce;
begin
setColor(black);
lineBlob(FxPos,FyPos,Length,Angle);  {Erase old}
if xForce=0 then xForce:=0.00001;
if yForce=0 then yForce:=0.00001;
angle:=arcTan(yForce/xForce);
if xForce>0 then angle:= angle+Pi;
FxPos:= round(beadX + beadR * cos(angle));
FyPos:= round(beadY + beadR * sin(angle));
length:=round(100 * sqrt(sqr(xForce)+sqr(yForce)));
setColor(red);
lineBlob( FxPos,FyPos,length,angle);
end;

procedure showAxes;
begin
setColor(blue);
line(0,yMid,xMax,yMid);  line(xMid,0,xMid,yMax);
end;

procedure showBead;
begin
setColor(black);
circle(xPos,yPos,beadR);  {erase old circle}
xPos:=round(beadX); yPos:=round(beadY);
setColor(cyan);
circle(xPos,yPos,beadR);  {draw new}
end;

procedure launchPhotonRight (ray:integer);
var  angle,xstart,ystart:real;  i:integer;
begin
  with photon[ray] do
  begin
  color:=yellow;
  angle:=spread*(random-0.5);  {random angles}
  Px:=1; Py:=0;  {set momenta for horizontal path toward right if colimated}
  if not colimated then begin Px :=cos(angle); Py := sin(angle); end;
  xStart:= xMid * (1-cos(angle));
  yStart:=yMid - xMid * sin(angle) + yOffset;
  for i:=0 to pLength do
    begin
    track[i,1] := xStart + Px * i;
    track[i,2] := yStart + Py * i;
    putPixel(round(track[i,1]),round(track[i,2]),color);
    end;
  endPtr:=pLength;
  inHiIndex:=false;
  end;{with photon[ray]}
end;

procedure launchPhotonLeft (ray:integer);
var  xStart, yStart, angle:real;  i:integer;
begin
  with photon[ray] do
  begin
  color:= green;
  angle:=Pi + spread * (random-0.5);
  Px:=-1; Py:=0;
  if not colimated then begin Px :=cos(angle);Py := sin(angle); end;
  xStart:= xMid * (1 - cos(angle));
  yStart:=yMid - xMid * sin(angle);
  for i:=0 to pLength do
    begin
    track[i,1] := xStart + Px * i;
    track[i,2] := yStart + Py * i;
    putPixel(round(track[i,1]),round(track[i,2]),color);
    end;
  endPtr:=pLength;
  inHiIndex:=false;
  end;{with photon[]}
end;

procedure limit (var arg:real; max,min:real);
begin
if arg>max then arg:=max;
if arg<min then arg:=min;
end;

procedure moveTheBead;  {apply external force to bead}
var correction:real;
begin  {assume strong overdamping and implicit time step per iteration}
targetX:=targetX - xforce;  {apply external force to bead}
targetY:=targetY - yforce;
{targetX, targetY also modified in "propagatePhotons" by light force}
correction:= targetX-beadX;
limit(correction,0.5,-0.5); {bead must not move faster than light}
beadX:=beadX+correction;
correction:= targetY-beadY;
limit(correction,0.5,-0.5);  {or else explosions occur}
beadY:=beadY+correction;
showBead;
end;

procedure propagatePhotons;
var oldHeadX,oldheadY,normX,normY,normL,oldPx,oldPy,
    normPx,normPy,tangPx,tangPy:real;
    i:integer;
    insideBead:boolean;
begin
for i:=-rayNum to rayNum do
    with photon[i] do
  begin
  oldHeadX:=track[endPtr,1]; oldHeadY:=track[endPtr,2];
  inc(endPtr); if endPtr>pLength then endPtr:=0; {move pointer to tail}
  putPixel(round(track[endPtr,1]), round(track[endPtr,2]), black); {erase it}
  track[endPtr,1]:= oldHeadX + Px; {make new head}
  track[endPtr,2]:= oldHeadY + Py;
  putPixel(round(track[endPtr,1]), round(track[endPtr,2]), color); {show it}
  insideBead := (sqr(track[endPtr,1]-beadX)+ sqr(track[endPtr,2]-beadY))
                 < sqr(beadR);
  if insideBead <> inHiIndex then  {photon hits bead/water interface}
      begin
      oldPx:=Px; oldPy:=Py;
      normL:=sqrt(sqr((track[endPtr,1]-beadX))+sqr(track[endPtr,2]-beadY));
      normX:=(track[endPtr,1]-beadX)/normL; {get unit normal to bead surface}
      normY:=(track[endPtr,2]-beadY)/normL;
      normPx:=normX * (normX*Px + normY*Py); {projection of P normal to surface}
      normPy:=normY * (normX*Px + normY*Py);
      tangPx:=Px-normPx; {projection of P tangential to surface}
      tangPy:=Py-normPy;
      inc(photoCount);
      if (photoCount mod transmitRate) = 0 then
         begin {photon is reflected}
         Px:=tangPx-normPx;
         Py:=tangPy-normPy;
         color:=white;  inc(reflCtr);
         end
      else
         begin {photon is refracted}   inc(refrCtr);
         inHiIndex:=insideBead; {adjust to new index}
         {tangPx and tangPy stay the same by Snell's law}
         if inHiIndex {but magnitude of P increases from 1 to refInd}
            then begin
            normPx:=normX * sqrt(sqr(relInd)-sqr(tangPx)-sqr(tangPy));
            normPy:=normY * sqrt(sqr(relInd)-sqr(tangPx)-sqr(tangPy));
            Px:=-normPx+tangPx;
            Py:=-normPy+tangPy;
            end;
         if (not inHiIndex) and ((sqr(tangPx)+sqr(tangPy))<1)
            then begin   {magnitude of P decreases to one}
            normPx:=normX * sqrt(1-sqr(tangPx)-sqr(tangPy));
            normPy:=normY * sqrt(1-sqr(tangPx)-sqr(tangPy));
            Px:= normPx+tangPx;
            Py:= normPy+tangPy;
            end;
         end;
      targetX:=targetX+LightUmph*(oldPx-Px);{apply light force to bead}
      targetY:=targetY+LightUmph*(oldPy-Py);{due to momentum change}
      end;{if photon hits interface}
  end;  {for photon[ray]}
end;

procedure showMenu;
var left, mid, right, height:integer;
begin
clearDevice;
left:=10; mid:=xMax div 3; right:=2*xMax div 3; height:=yMax-60;
setColor(lightBlue);
outTextXY(left,height,'m_____menuOff');
outTextXY(mid,height,'arrows__applyForces');
outTextXY(right, height,'s__singleBeam');
height := height + 15;
outTextXY(left,height,'w_____widenBeams');
outTextXY(mid, height,'n_______narrowBeams');
outTextXY(right,height,'d__dualBeams');
height:=height+15;
outTextXY(left,height,'home__raiseBeam');
outTextXY(mid,height,'end_____lowerBeam');
outTextXY(right,height,'c__colimateBeam');
height:=height+15;
outTextXY(left,height,'h____haltPhotons');
outTextXY(mid,height,'z_____zeroBead');
outTextXY(right,height,'q__quit');
end;

procedure startScreen;
begin
clearDevice;
showBead;
setColor(lightBlue);
outTextXY(10,yMax-40,'m__menu');
outTextXY(10,yMax-30,'q__quit');
end;

procedure zeroBead;
var i:integer;
begin
beadX:=xMid; beadY:=yMid;
targetX:=beadX; targetY:=beadY;
for i:=-rayNum to rayNum do with photon[i] do
    inHiIndex:= (sqr(track[endPtr,1]-beadX) +sqr(track[endPtr,2]-beadY))
     < sqr(beadR);
startscreen;
end;

begin {main}
grDriver:=detect;
initgraph(grDriver,grMode,'');
xMax:=getMaxX; yMax:=getMaxY; xMid:=xMax div 2; yMid:=yMax div 2;
beadX:=xMid; beadY:=yMid;
targetX:=beadX; targetY:=beadY;
beadR:=50; spread:=1.0; colimated:=false;
xforce:=0.0;yforce:=0.0;
transmitRate:=round(1/sqr((relInd-1)/(relInd+1))); {according to Fresnel
eqn., transmitRate = number of photons that will transmit before one refelcts.}
startScreen;
repeat
if (loopCtr mod 4 = 0) then launchPhotonRight(loopCtr div 4);
if (loopCtr mod 4 = 2) then if dualBeam
            then launchPhotonLeft(-(loopCtr div 4))
            else launchPhotonRight(-(loopCtr div 4));
propagatePhotons;
moveTheBead;
if loopCtr mod 10 = 1 then begin showBead; showAxes; showForce; end;
inc(loopCtr);
if loopCtr>800 then loopCtr:=0;
if (loopCtr mod 20 = 1) and keyPressed then
  begin
  ch:=readkey;
     case ch of
     's':dualBeam:=false;
     'd':dualBeam:=true;
     'w':spread:=spread*1.1;
     'n':spread:=spread/1.1;
     'c':colimated:=not colimated;
     'r':begin  outTextXY(xMid+5,20,'what fractional reflection?');
         readln(rate); rate:=rate+1e-4; {prevent div by zero}
         transmitRate:=round(1/rate);clearDevice;
         end;
     'z':zeroBead;
     'h':ch:=readkey;
     'b':clearDevice;
     'm': begin menuOn := not menuOn;
          if menuOn then showMenu else startScreen;
          end;
      chr(0): begin ch:=readkey;
                case ord(ch) of
    {up arrow}  72: yforce:=yforce+0.05;
{down arrow}    80: yforce:=yforce-0.05;
{right arrow}   77: xforce:=xforce-0.05;
 {left arrow}   75: xforce:=xforce+0.05;
 {home}         71: dec(yOffset);
 {end}          79: inc(yOffset);
                end; {inner case}
              end; {chr(0)}
     end;{outer case}
  end; {keypressed}
until (ch = 'q');
end.