clear;

format long;

DMAX = 10000;

VMAX = 0.01;

N = 6;

alpha = 7.28036e-3;

theta = 1.45612e-2;

K = (1.602176e-19^2) / (4*pi*8.854e-12*9.109382e-31*2.99792458e8^2);

dv = VMAX/DMAX;

v = dv;

te = asin(v);

chkcross = 1;

r = 1;

for i = 1:DMAX;

solution(i) = (sin((v/alpha)*((N*pi)+((-1^N)*theta)))/v) - 1;

radii(i) = K / (sin(te) * tan(te));

vel(i) = v;

if chkcross == 1

if solution(i) < 0

permissible_orbits(r, 1) = radii(i);

permissible_orbits(r, 2) = vel(i);

r = r+1;

chkcross = -1;

end

end

if chkcross == -1

if solution(i) > 0

permissible_orbits(r, 1) = radii(i);

permissible_orbits(r, 2) = vel(i);

r = r+1;

chkcross = 1;

end

end

v = v+dv;

te = asin(v);

end

permissible_orbits

plot(vel, solution);

axis([0, 0.008, -600, 600]);