function null
% Schwarzchild null geodesic
% of the first and the second kind
P = 1;     
M = 2/14;  % mass of the black hole

Q = sqrt((P-2*M)*(P+6*M));
k = sqrt((Q-P+6*M)/2/Q);


chispan = [2*asin(sqrt((Q-P+2*M)/(Q-P+6*M))) 3.75];
                        %2*asin(sqrt((Q-P+2*M)/(Q-P+6*M)))
options = odeset('Refine',38);
[fi1, chi1] = ode45(@fi, chispan, 2.2, options, Q, ...
    P, k); 

chispan = [0 pi];
[fi2, chi2] = ode45(@fi, chispan, 0, options, Q, ...
    P, k);

u1 = 1/P -(Q-P+6*M)/8/M/P*(1+cos(chi1));
u2 = 1/P +(Q+P-6*M)/4/M/P*sec(chi2/2).^2;

polar(fi1, 1./u1)
hold on
%polar(-fi1, 1./u1)
polar(fi2, 1./u2)
title(['M = ' num2str(M) ', P = ' num2str(P)], 'Fontsize', ...
    15);
%----------------------------------------------

function f = fi(fi, chi, Q, P, k)
% ODE for phi coordinate for the null geodesics 
% f1 is d chi / d phi
f = sqrt(Q/P) *sqrt(1-k^2*sin(chi/2)^2);