(* Copyright (C)2024 Sasan Ardalan                                      *)
(* This program is free software:you can redistribute it and/or modify  *)(* it under the terms of the GNU General Public License as published by *)
(* the Free Software Foundation,either version 3 of the License,or      *)(*(at your option) any later version.This program is distributed        *)
(* in the hope that it will be useful,but WITHOUT ANY WARRANTY          *)
(* without even the implied warranty of                                 *)
(* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.See the          *)
(* GNU General Public License for more details.You should have          *)
(* received a copy of the GNU General Public License                    *)
(* along with this program.If not,see<http://www.gnu.org/licenses/>.   *)


(* Velocity Vector for Elliptic Orbits           *)
(* Orbital Mechanics                             *)
(* Author: Sasan Ardalan                         *)
(* ScientificWorks.org                           *)
(* Date: September 10, 2024                      *)



p=1.0
e=0.6
f[theta_]:=p/(1+e*Cos[theta])

PolarPlot[f[theta],{theta,0,2*Pi}]

dyDivdx[theta_]:=( e*(Sin[theta])^2+f[theta]*Cos[theta]*(1+e*Cos[theta])^2)/(e*Sin[theta]*Cos[theta]-f[theta]*Sin[theta]*(1+e*Cos[theta])^2)

(* PolarPlot[dyDivdx[theta],{theta,0,Pi/2}] *)

arrows={}
delTheta=0.1*Pi
theta=0;
For[i=0,i<20, i++,
theta=theta+delTheta;
(*theta=Pi*(0.1); *)

(* Print["Theta:",theta]; *)
m=dyDivdx[theta];
(* Print["Slope m:",m]; *)
r=f[theta];
x0=r*Cos[theta];
y0=r*Sin[theta];
(* Print["x0=",x0," y0=", y0]; *)
n=y0-m*x0;
ftan[x_]:=m*x+n;
x=x0;
yy=ftan[x0];
(* Print["x0==",x0," yy=", y0]; *)
(*Plot[ftan[x],{x,-4,2}]*)

lenF=0.4;
ra=1/(1+e);
rb=1/(1-e);
a=(ra+rb)/2;
(* Print["a first=",a]*)
a=1/(1-e^2);
(* Print["a=",a]; *)

velocity=Sqrt[(2*a-r)/(r*a)];
(*Print ["Velocity =",velocity]; *)

deltax=lenF/Sqrt[1+m^2]*velocity;
x11=If[m<0,x0-deltax,x0];
x22=If[m<0, x0, x0+deltax];

x1=x11;
x2=x22;
(*Print["x0=",x0,"   x1=",x1,"   x2=",x2];*)

arrows=If[theta > Pi, Append[arrows,Arrow[{{x1,ftan[x1]},{x2,ftan[x2]} }]],Append[arrows,Arrow[{{x2,ftan[x2]},{x1,ftan[x1]} }]]]


] 

Show[PolarPlot[f[theta],{theta,0,2*Pi}],Graphics[arrows]]
(*Print[arrows]*)