rem Ellipsenkaustik
rem Was machen die Strahlen, die von einem leuchtenden Punkt
rem im Innern einer verspiegelten Ellipse ausgehen? 

rem Bildgršsse
horz0=0: rem links 
vert0=0: rem oben
horz1=900: rem rechts
vert1=720: rem unten
rem printable pixels: horz0..horz1-1, vert0..vert1-1
black=0:white=65535
forecolor white,white,white: rem Vordergrund (pixelfarbe) weiss (= rot+grŸn+blau)
backcolor black,black,black: rem Hintergrund 
resize console horz0,vert0,horz1,vert1
virtu=init screen(horz0,vert0,horz1,vert1): rem virtueller Bildschirm 
set screen to virtu
cls
randomize timer
system breathe

pi=3.141592653589793238462643
zweipi=6.283185307179586476925287: rem 2*pi
pihalbe=1.570796326794896619231322: rem pi/2
def fn min(x,y)=-x*(x<=y)-y*(x>y)
def fn max(x,y)=-x*(x>=y)-y*(x<y)


rem Ellipsenparameter
a=(horz1-horz0)/2.1
b=(vert1-vert0)/2.1
c=sqr(a*a-b*b)
mx=(horz1+horz0)/2
my=(vert1+vert0)/2


rem Anzahl Strahlen, die von der Quelle ausgesandt werden
Anz_Strahlen=10000

rem Anzahl interne Reflexionen
anz_refl=2

rem Quellpunkt
xq=-1.1*c
yq=0.25*b
fcircle mx+xq,my+yq,2: rem ausgefŸllter Kreis mit Radius 2 Pixel

textsize 24
text 0,24," n = "+str$(anz_Strahlen)
text 0,48," r = "+str$(anz_refl)
text 0,72," x = 1.1 c"
text 0,96," y = 0.25 b"

rem Ellipse zeichnen
for x=0 to a step 1
   y=sqr(1-x*x/a/a)*b
   fcircle mx+x,my+y,1
   fcircle mx+x,my-y,1
   fcircle mx-x,my+y,1
   fcircle mx-x,my-y,1
next x
for y=0 to b step 1
   x=sqr(1-y*y/b/b)*a
   fcircle mx+x,my+y,1
   fcircle mx+x,my-y,1
   fcircle mx-x,my+y,1
   fcircle mx-x,my-y,1
next y

rem Strahlen zeichnen
for phi=0 to zweipi step zweipi/Anz_Strahlen 
   kx=cos(phi): rem Richtungsvektor
   ky=sin(phi)
   al=b*b*kx*kx+a*a*ky*ky
   be=b*b*2*xq*kx+a*a*2*yq*ky
   ga=b*b*xq*xq+a*a*yq*yq-a*a*b*b
   t=(-be+sqr(be*be-4*al*ga))/2/al
   xp=xq+kx*t: rem Auftreffpunkt auf der Ellipse
   yp=yq+ky*t
   for tt=0 to t step 1
      x=xq+kx*tt: rem vom Quellpunkt zum Punkt auf der Ellipse
      y=yq+ky*tt
      get pixel mx+x,my+y,red,green,blue
      red=fn min(red+256,65535): rem Grauwert des Pixels erhšhen
      forecolor red,red,red
      plot mx+x,my+y
   next tt
   rem reflektierte Strahlen zeichnen
   for i=1 to anz_refl
      sx=-xp/a/a: rem Normale im Auftreffpunkt
      sy=-yp/b/b
      s=sqr(sx*sx+sy*sy)
      sx=sx/s
      sy=sy/s
      ks=kx*sx+ky*sy
      kx=kx-sx*ks*2: rem Richtungsvektor reflektierter Strahl
      ky=ky-sy*ks*2
      al=b*b*kx*kx+a*a*ky*ky
      be=b*b*2*xp*kx+a*a*2*yp*ky
      ga=b*b*xp*xp+a*a*yp*yp-a*a*b*b
      t=(-be+sqr(be*be-4*al*ga))/2/al
      xr=xp+kx*t: rem Auftreffpunkt des reflektierten Strahls auf die Ellipse
      yr=yp+ky*t
      for tt=0 to t step 1
         x=xp+kx*tt: rem reflektierten Strahl zeichnen
         y=yp+ky*tt
         get pixel mx+x,my+y,red,green,blue
         red=fn min(red+256,65535)
         forecolor red,red,red
         plot mx+x,my+y
      next tt
      xp=xr
      yp=yr
   next i
   if button then goto ausgang: rem Abbruch bei Mausklick
next phi


ausgang: 

forecolor white,white,white
set screen to console
rem virtuellen Bildschirm auf Screen kopieren (zum anzeigen und speichern)
copyrect horz0,vert0,horz1,vert1,horz0,vert0,horz1,vert1,0,virtu,0
rem Screen als .pict-Bild speichern
save pict "Ellipsenkaustik p07-05-1",horz0,vert0,horz1,vert1

print "fertig"
kill screen virtu
end