%%  Lichtablenkung

% W. Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB 
% 2016

figure
r0=5;                           % ist nicht der minimale Abstand
rs = 1;                         % rs Schwarzschildradius
rmin = 1/(1/r0+rs/(2*r0^2));
phi0 = -90.0/180*pi;            % Startwinkel
                                % Anfangswerte zur Loesung der DGL aus
                                % Naeherungsrechnung
y00 = cos(phi0)/r0 + rs/(2*r0^2) * (1 + sin(phi0)^2);
dy0 = -sin(phi0)/r0 + rs/(r0^2) * sin(phi0)*cos(phi0);
y0 =[y00;dy0];

[phi,yy] = ode45(@(phi,y) LichtablenkungDGL(phi,y,rs), [phi0,pi/2], y0);  % DGL

% Visualisierung der Lichtablenkung
r = 1./yy(:,1);
x = r.*sin(phi);
y = r.*cos(phi);

plot(x,y),axis equal,shg
hold on
plot(x,-y)

% Visulaisierung der scheinbaren Sternorte
theta = linspace(0,2*pi);
area(rmin/2*sin(theta),rmin/2*cos(theta))%, grid on
yseh = (y(end)-y(end-1))/(x(end)-x(end-1)) *x - (y(end)-y(end-1))/(x(end)-x(end-1))*x(end);
plot(x,yseh,x(1),yseh(1),'p')
yseh = (y(end)+y(end-1))/(x(end)-x(end-1)) *x - (y(end)+y(end-1))/(x(end)-x(end-1))*x(end);
plot(x,yseh,x(1),yseh(1),'p')
plot(x(1),y(1),'o')
hold off
%hold on
%plot([0,0],[0,r0],'o')
