//---------------------------------------------------------------------------
//	Projekcja relatywistyczna podstawowych form geometrycznych

#include <vcl\vcl.h>
#pragma hdrstop

#include "natura.h"				//prędkość światła w układzie SI
#include "relatyw3d.h"
#include "math.h"
//---------------------------------------------------------------------------
// Konstruktor główny, inicjujący wszystkie parametry przekształcenia 3d
TRelatyw3d :: TRelatyw3d(
   TPunkt Aobserwator,                     //pozycja obserwatora
   int xe0, int ye0, int eszer, int ewys,  //okno ekranowe
   double v,                               //prędkość w kierunku osi 'x'
   TCanvas *Acnv,                          //aparat kreślący
   double Aodl_ekr/*=0.5*/,
   double szer_ekr/*=0*/, double wys_ekr/*=0*/)//fizyczny opis ekranu

{
 double r1, a;

 if( szer_ekr <=0)                          //gdy nie określono rozmiarów ekranu
    szer_ekr = (double)eszer * 0.0257 / 96.;//0.0257 m = 1 cal. Ekran ma 96 pikseli/cal
 if( wys_ekr <=0)
    wys_ekr  = (double)ewys  * 0.0257 / 96.;


 beta = v / N_c;                            //parametry kinematyki relatywistycznej
 if( fabs( beta) < 1.0)
 	gamma = 1. / sqrt( 1. - beta * beta);
 else
 	gamma = 1e10;                           //osiągnięto prędkość światła ...


 obserwator = Aobserwator;                  //przechowanie niezbędnych parametrów
 odl_ekr = Aodl_ekr;
 cnv = Acnv;


 a = obserwator.x * obserwator.x + obserwator.y * obserwator.y;
 r1 = sqrt( a);                             //odl. w płaszczyźnie XY
 odl_obserwatora = sqrt( a + obserwator.z * obserwator.z);//odl. obserwatora od centrum
 if( r1 != 0)                               //wyliczenie współczynników macierzy obrotu
 {
    sf = obserwator.y / r1;                 //sin( fi)
    cf = obserwator.x / r1;                 //cos( fi)
 }
 else
 {
    sf = 0;
    cf = 1;
 }
 st = r1 / odl_obserwatora;                 //sin( teta)
 ct = obserwator.z / odl_obserwatora;       //cos( teta)

 A = (double)eszer / szer_ekr;              //przekszt. lin.
 B = xe0 + A * szer_ekr / 2.;
 C = -(double)ewys / wys_ekr;
 D = ye0 - C * wys_ekr / 2.;
}
//---------------------------------------------------------------------------
//  Wylicz współrzędne punktu relatywistycznego.
//  Jest to najważniejszy algorytm niniejszego zagadnienia.
//  Tutaj znajdują się receptury przekształcenia relatywistycznego.
bool TRelatyw3d :: wylicz_punkt_relatyw( int &xe, int &ye, TPunkt p)
{
 double xn, yn, zn;                     //punkt 3d w układzie komputerowca
 double odl_punktu;                     //odległość punktu od obserwatora
 const int max = 10000;                 //Grafika się niekiedy wiesza ...

 odl_punktu = sqrt(( p.x - obserwator.x) * ( p.x - obserwator.x) +
                   ( p.y - obserwator.y) * ( p.y - obserwator.y) +
                   ( p.z - obserwator.z) * ( p.z - obserwator.z));

//-------------------------------------
 p.x = p.x / gamma;                     //skrócenie relatywistyczne ...
 p.x = p.x - ( odl_punktu - odl_obserwatora) * beta;// ...i opóźnienie sygnału, odniesione do środka układu
//-------------------------------------

 xn = - 	 sf * p.x + cf * p.y;       //transformacja
 yn = - cf * ct * p.x - sf * ct * p.y + st * p.z;
 zn = - cf * st * p.x - sf * st * p.y - ct * p.z + odl_obserwatora;

 if( zn == 0)                           //punkt drapnął komputerowca w oko ...
	return false;

 xe = A * xn * odl_ekr / zn + B;        //twierdzenie Talesa ...
 if( xe < -max)
    xe = -max;
 if( xe > max)
    xe = max;

 ye = C * yn * odl_ekr / zn + D;
 if( ye < -max)
    ye = -max;
 if( ye > max)
    ye = max;

 return true;
}
//---------------------------------------------------------------------------
//  Wykreśl punkt relatywistyczny
bool TRelatyw3d :: punkt_relatyw( TPunkt p, TColor kolor/*=clBlack*/)
{
 int xe, ye;
 bool ret;

 ret = wylicz_punkt_relatyw( xe, ye, p);

 if( ret)
 	cnv -> Pixels[ xe][ ye] = kolor;

 return ret;
}
//---------------------------------------------------------------------------
//  Wykreśl linię relatywistyczną.
//  Uwaga: linia może mieć obraz powyginany - trzeba ją kreślić małymi odcinkami!
//  Każdy odcinek może mieć inny kolor
bool TRelatyw3d :: linia_relatyw( TPunkt p1, TPunkt p2, int il_podzial,
    void (__closure* kolor_fun)( TCanvas *, int)/*=NULL*/)
{
 int i, xe, ye;
 TPunkt p;
 TLinia l(p1, p2);
 bool ret = true;

 for( i = 0; i <= il_podzial; i ++)
 {
    p = l.wylicz_punkt_na_linii( i, il_podzial);
    if( wylicz_punkt_relatyw( xe, ye, p))
    {
        if( i == 0)
            cnv -> MoveTo( xe, ye);
        else
        {
            if( kolor_fun)                  //jeśli określono funkcję koloru ...
                kolor_fun( cnv, i);
            cnv -> LineTo( xe, ye);
        }
    }
    else
    {
        ret = false;                        //nie udało się dokreślić linii do końca ...
        break;
    }
 }
 return ret;
}
//---------------------------------------------------------------------------
//  Wykreśl linię relatywistyczną.
//  Uwaga: linia może mieć obraz powyginany - trzeba ją kreślić małymi odcinkami!
//  Cała linia ma ten sam kolor
bool TRelatyw3d :: uproszczona_linia_relatyw( TPunkt p1, TPunkt p2,
                int il_podzial, TColor kolor)
{
 int i, xe, ye;
 TPunkt p;
 TLinia l(p1, p2);
 bool ret = true;

 cnv -> Pen -> Color = kolor;
 for( i = 0; i <= il_podzial; i ++)
 {
    p = l.wylicz_punkt_na_linii( i, il_podzial);
    if( wylicz_punkt_relatyw( xe, ye, p))
    {
        if( i == 0)
            cnv -> MoveTo( xe, ye);
        else
            cnv -> LineTo( xe, ye);
    }
    else
    {
        ret = false;                        //nie udało się dokreślić linii do końca ...
        break;
    }
 }
 return ret;
}
//---------------------------------------------------------------------------
//  Wykreśl czworokąt relatywistyczny.
//  Uwaga: płaski czworokąt może mieć obraz powyginany - trzeba  kreślić kawałkami!
//  Przed wykreśleniem każdego elementu siatki wywoływana jest funkcja
//  użytkownika kolor_fun( int i, int j), pozwalająca ustalić atrybuty
//  (np. kolor) elementu (i, j).
bool TRelatyw3d :: czworokat_relatyw( TPunkt p1, TPunkt p2, TPunkt p3, TPunkt p4,
        int il_podzial,
        void (__closure* kolor_fun)( TCanvas *, int, int)/*=NULL*/)
{
 int i, j, xe, ye;
 TPunkt **wezly;                            //węzły siatki
 TPoint  **ekran;                            //te same węzły, ale na ekranie
 TPoint p[ 4];                               //elementarny czworokącik nakładanej siatki
 TLinia l1( p1, p4), l2( p2, p3);           //brzegi obszaru (boki drabiny)
 TLinia l;                                  //linia rozpinana m-dzy brzegami obszaru

 if( il_podzial <= 0)
    il_podzial = 1;

 try
 {
    wezly = new TPunkt* [ il_podzial + 1];  //tablica węzłów siatki, nakładanej na czworokąt
    ekran = new TPoint * [ il_podzial + 1];  //... i ich reprezentacji ekranowych
    for( i = 0; i <= il_podzial; i ++)
    {
        wezly[ i] = new TPunkt [ il_podzial + 1];
        ekran[ i] = new TPoint  [ il_podzial + 1];
    }
 }
 catch(...)
 {
    return false;
 }

 for( i = 0; i <= il_podzial; i ++)         //wylicz punkty węzłowe
 {
    l = TLinia( l1.wylicz_punkt_na_linii( i, il_podzial),//szczeble drabiny
                l2.wylicz_punkt_na_linii( i, il_podzial));

    for( j = 0; j <= il_podzial; j ++)
    {
        wezly[ i][ j] = l.wylicz_punkt_na_linii( j, il_podzial);
    }
 }


 for( i = 0; i <= il_podzial; i ++)     //wylicz współrzędne ekranowe węzłów
 {
    for( j = 0; j <= il_podzial; j ++)
    {
        wylicz_punkt_relatyw( xe, ye, wezly[ i][ j]);
        ekran[ i][ j].x = xe;
        ekran[ i][ j].y = ye;
    }
 }

 for( i = 0; i < il_podzial; i ++)      //wykreśl powierzchnię
 {
    for( j = 0; j < il_podzial; j ++)
    {
        p[ 0] = ekran[ i][ j];
        p[ 1] = ekran[ i][ j + 1];
        p[ 2] = ekran[ i + 1][ j + 1];
        p[ 3] = ekran[ i + 1][ j];
        if( kolor_fun)                  //jeśli określono funkcję koloru ...
            kolor_fun( cnv, i, j);
        cnv -> Polygon( p, 3);
    }
 }


 for( i = 0; i <= il_podzial;  i++)     //likwidacja już niepotrzebnych współrzędnych ekranowych
 {
    delete[] ekran[ i];
 }
 delete [] ekran;

 for( i = 0; i <= il_podzial;  i++)     //likwidacja już niepotrzebnej tablicy węzłów
    delete[] wezly[ i];
 delete [] wezly;

 return true;
}
//---------------------------------------------------------------------------
//  Wykreśl czworokąt relatywistyczny.
//  Uwaga: płaski czworokąt może mieć obraz powyginany - trzeba  kreślić kawałkami!
//  Wszystkie kawałki mają ten sam kolor.
bool TRelatyw3d :: uproszczony_czworokat_relatyw( TPunkt p1, TPunkt p2, TPunkt p3, TPunkt p4,
        int il_podzial, TColor k1, TColor k2)
{
 int i, j, xe, ye;
 TPunkt **wezly;                            //węzły siatki
 TPoint  **ekran;                            //te same węzły, ale na ekranie
 TPoint p[ 4];                               //elementarny czworokącik nakładanej siatki
 TLinia l1( p1, p4), l2( p2, p3);           //brzegi obszaru (boki drabiny)
 TLinia l;                                  //linia rozpinana m-dzy brzegami obszaru

 if( il_podzial <= 0)
    il_podzial = 1;

 try
 {
    wezly = new TPunkt* [ il_podzial + 1];  //tablica węzłów siatki, nakładanej na czworokąt
    ekran = new TPoint * [ il_podzial + 1];  //... i ich reprezentacji ekranowych
    for( i = 0; i <= il_podzial; i ++)
    {
        wezly[ i] = new TPunkt [ il_podzial + 1];
        ekran[ i] = new TPoint [ il_podzial + 1];
    }
 }
 catch(...)
 {
    return false;
 }

 cnv -> Pen -> Color = k1;
 cnv -> Brush -> Color = k1;
 for( i = 0; i <= il_podzial; i ++)         //wylicz punkty węzłowe
 {
    l = TLinia( l1.wylicz_punkt_na_linii( i, il_podzial),//szczeble drabiny
                l2.wylicz_punkt_na_linii( i, il_podzial));

    for( j = 0; j <= il_podzial; j ++)
    {
        wezly[ i][ j] = l.wylicz_punkt_na_linii( j, il_podzial);
    }
 }


 for( i = 0; i <= il_podzial; i ++)     //wylicz współrzędne ekranowe węzłów
 {
    for( j = 0; j <= il_podzial; j ++)
    {
        wylicz_punkt_relatyw( xe, ye, wezly[ i][ j]);
        ekran[ i][ j].x = xe;
        ekran[ i][ j].y = ye;
    }
 }

 for( i = 0; i < il_podzial; i ++)      //wykreśl powierzchnię
 {
    for( j = 0; j < il_podzial; j ++)
    {
        p[ 0] = ekran[ i][ j];
        p[ 1] = ekran[ i][ j + 1];
        p[ 2] = ekran[ i + 1][ j + 1];
        p[ 3] = ekran[ i + 1][ j];
        cnv -> Polygon( p, 3);
    }
 }


 for( i = 0; i <= il_podzial;  i++)     //likwidacja już niepotrzebnych współrzędnych ekranowych
 {
    delete[] ekran[ i];
 }
 delete [] ekran;

 for( i = 0; i <= il_podzial;  i++)     //likwidacja już niepotrzebnej tablicy węzłów
    delete[] wezly[ i];
 delete [] wezly;

 return true;
}

