/*
 Kaustik mit der Lieberherr-Methode
 in geschlossenen Kegelschnitten (kreis,Ellipse)
 
 Autor
 Daniel Baumgartner
 Kantonsschule Baden
 
 www.mathematicus.ch
 eMail <vorname>_<nachname>@bluewin.ch
*/ 
 
#include <iostream>
#include <cmath>
using namespace std;

#include "gsm.h"

/*Die Koeffizienten des Kegelschinitts
  k: k1*x^2+k2*y^2+k3*z^2+k4*y*z+k5*x*z+k6*x*y=0 (z:=1)
*/

/* Parameter der Ellipse
 a, b die Laengen der Halbachsen
*/
double a=2.0;
double b=2*0.681;

double k1=b*b;
double k2=a*a;
double k3=-a*a*b*b;
double k4=0;
double k5=0;
double k6=0;

// Brennpunkt der Ellipse
double f=sqrt(a*a-b*b);

/* 
  Die Lichtquelle
*/

// pkt L(2.0,0.0);
// pkt L(f,0.0);
pkt L(a/2,b/2);

/*
 Berechnung des neben dem Reflexionspunkt zweiten Schnittpunkts 
 des reflektierten Strahls
*/

pkt reflexions_punkt(double t)
{
    pkt h;
    h.def(a*cos(t),b*sin(t));
    return h;
};

vek tangenten_vektor(double t)
{
    vek v;
    v.def(-a*sin(t),b*cos(t));
    return v;
};

pkt end_punkt(pkt L,double t)
{
    pkt e;
    pkt p=reflexions_punkt(t);
    vek u=tangenten_vektor(t);
    // Lot f={f1,f2,f3} in P
    double f1=u.x;
    double f2=u.y;
    double f3=-p.x*u.x-p.y*u.y;
    // L*={l1,l2,l3} der an f gespiegelte Punkt L
    double l1=2*f1*f3+f1*f1*L.x-f2*f2*L.x+2*f1*f2*L.y;
    double l2=2*f2*f3+2*f1*f2*L.x-f1*f1*L.y+ f2*f2*L.y;
    double l3=-f1*f1-f2*f2;
    // die Gerdade PL*
    double g1=-l2+l3*p.y;
    double g2=l1-l3*p.x;
    double g3=l2*p.x-l1*p.y;
    // der neben P zweite Schnittpunkt der Geraden PL* mit dem Kegelschnitt
    double p1=-2*g1*g3*k2+g1*g2*k4-g2*g2*k5+g2*g3*k6-g2*g2*k1*p.x-g1*g1*k2*p.x+g1*g2*k6*p.x;
    double p2=-2*g2*g3*k1-g1*g1*k4+g1*g2*k5+g1*g3*k6-g2*g2*k1*p.y-g1*g1*k2*p.y+g1*g2*k6*p.y;
    double p3=g2*g2*k1+g1*g1*k2-g1*g2*k6;
    e.def(p1/p3,p2/p3);
    return e;
};

/*
 Hauptteil
*/
 
main()
{
// Koordinatenbereich xmin,xmax,ymin,ymax
pixel_koord(-2.2,2.2,-2.2,2.2);

pix_X=1024;//Anzahl Pixel in x-Richtung

// Die Matrix festlegen
int m=(int) round(deltay/deltax*pix_X);
intmat pix(pix_X,m);

    
int anz=1024;//Anzahl Lichtstrahlen
double w=2*pi/(double)anz;
pkt p;

for(long i=1;i<anz;i++){
    p=end_punkt(L,w*i);
    pixel_strecke(L,reflexions_punkt(w*i),pix);
    pixel_strecke(reflexions_punkt(w*i),p,pix);
};

/* Die Graustufenmatrix erzeugen
 0=Mathematica-Syntax
 1=Octave-/Matlab-Syntax
*/
pixelbild(pix,"ellipse02",1);

}