%variable gravity in a sphere
Rc = 471000;%radius of ceres
Rm = 1737100;%radius of moon
R = 3389500;%radius of mars
Re = 6371000;%radius of earth
G = 6.67384*10^(-11);%gravitational constant
Mc = 9.43*10^20;%mass of ceres
Mm = 7.3477*10^22;%mass of moon
M = 6.4174*10^23;%mass of mars
Me =4.9736*10^24;%mass of earth
a = 6.022141*10^23;%avagadros number
%m represents kg/particle of the atmosphere
%m = 1*10^-3/a;%mass of hydrogen kg/particle
%m = 18.01528*10^-3/a;% mass of water
%m = 15.9994*10^-3/a;% mass of O2
m = 43.34*10^-3/a;%average molecular weight currently on mars kg/particle
me = 28.97*10^-3/a;% average molecular weight on earth kg/particle
%me = 28*10^-3/a;% molar mass of N2
%n0 represents particle density at the surface of a planet
n0e = 1.217/me;%1.217 kg/m^3 is the density on earth's surface
n0 = .02/m;%.02 kg/m^3 is the density on mars's surface
T = 210;%average tempurature on mars now
Te = 288;%average tempurature on earth
k = 1.3806488*10^-23;%boltzmann constant
B = 1/(k*T);%beta
Be = 1/(k*Te);
y1 = linspace(0,111000, 10000);
n = n0*exp(B*G*M*m./(y1+R))./exp(B*G*m*M/R);%particle density array
%ne = n0e*exp(Be*G*Me*me./(y1+Re))./exp(Be*G*me*Me/Re);%particle density array for earth
nme =  n0e*exp(Be*G*M*me./(y1+R))./exp(Be*G*me*M/R);%particle density array for mars with earthlike conditions
nmoe =  n0e*exp(Be*G*Mm*me./(y1+Rm))./exp(Be*G*me*Mm/Rm);%particle density array for mars with earthlike conditions
nve =  n0e*exp(Be*G*Mc*me./(y1+Rc))./exp(Be*G*me*Mc/Rc);%particle density array for mars with earthlike conditions
figure(1)% density model
hold off
subplot(2,1,1)
semilogx(n*m,y1,'g')
axis([10e-16 1 0 1000000])
ylabel('altitude (m)');
xlabel('atmospheric density at surface kg/m^3');
title('semilog plot of the atmospheric model of Mars');

subplot(2,1,2)
plot(n*m,y1,'g')

hold on
%plot(ne*me,y1,'b')
%plot(nme*me,y1,'r')
%plot(nmoe*me,y1,'y')
%plot(nve*me,y1,'b')
%legend('current Mars model','Mars with Earth-like atmosphere')
ylabel('altitude (m)');
xlabel('atmospheric density at surface kg/m^3');
title('atmospheric model of Mars');
Vt = sqrt(3*k*T/(m))%average thermal velocity
Vtn = sqrt(3*k*Te/(me))%average thermal velocity earthlike conditions
y2 = linspace(0,1000000000,10000);
Ve = sqrt(2*G*M./(R+y2));%escape velocity array based on y
Ven = sqrt(2*G*M./(R+y2));%escape velocity array based on y earthlike conditions
Vem = sqrt(2*G*Mm./(Rm+y2));%escape velocity array based on y earthlike conditions moon
Vev = sqrt(2*G*Mc./(Rc+y2));%escape velocity array based on y earthlike conditions ceres
Hc = 2*m*M*G/(3*k*T)-R%height at which Vt = Ve (m)
Hcn = 2*me*M*G/(3*k*Te)-R%height at which Vt = Ve (m) earthlike conditions
Hcm = 2*me*Mm*G/(3*k*Te)-Rm%height at which Vt = Ve (m) earthlike conditions moon
Hcv = 2*me*Mc*G/(3*k*Te)-Rc%height at which Vt = Ve (m) earthlike conditions ceres
nHc = m*n0*exp(B*G*M*m./(Hc+R))./exp(B*G*m*M/R)%density at Hc kg/^3
nHcn = me*n0*exp(Be*G*M*me./(Hcn+R))./exp(Be*G*me*M/R)%density at Hc kg/^3 earthlike conditions
nHcm = me*n0*exp(Be*G*Mm*me./(Hcm+Rm))./exp(Be*G*me*Mm/Rm)%density at Hc kg/^3 earthlike conditions moon
nHcv = me*n0*exp(Be*G*Mc*me./(Hcv+Rc))./exp(Be*G*me*Mc/Rc)%density at Hc kg/^3 earthlike conditions ceres
figure (2)%escape velocity
plot(Ve,y2)
xlabel('escape velocity (m/s)');
ylabel('altitude (m)');
title('escape velocity on mars');
Vj = linspace(1,10000,10000);
Pdv = 4*pi*(m/(2*pi*k*T))^(3/2).*Vj.^2.*exp(-m.*Vj.^2/(2*k*T));% differential probability of having velocity Ve
Pdvn = 4*pi*(me/(2*pi*k*Te))^(3/2).*Vj.^2.*exp(-me.*Vj.^2/(2*k*Te));% differential probability of having velocity Ve earthlike conditions
figure (3)%differential probability
plot (Vj, Pdv)
title('differential probability of particle velocity given T = 210 K');
xlabel('velocity (m/s)');
ylabel('probability');
Pv = fliplr(cumsum(fliplr(Pdv)));%probability of finding a particle with velocity greater than v
Pvn = fliplr(cumsum(fliplr(Pdvn)));%probability of finding a particle with velocity greater than v earthlike conditions
figure (4)%integrated probability
plot (Vj, Pv)
title('probability of velocity greater than v given T = 210 K');
xlabel('velocity (m/s)');
ylabel('probability');
y3 = 2*G*M./Vj.^2-R;%escape height array for given velocity
y3m = 2*G*Mm./Vj.^2-Rm;%escape height array for given velocity moon
y3v = 2*G*Mc./Vj.^2-Rc;%escape height array for given velocity ceres
figure (5)%probability of escape
%plot(Pv, y3)
hold on
plot(Pvn, y3,'r')
plot (Pvn, y3m, 'y')
plot (Pvn, y3v, 'b')
hold off
legend('mars','moon', 'ceres')
axis([0 1 0 1000000000])
title('probability of a particle reaching escape velocity');
xlabel('probability');
ylabel('altitude (m)');
np =  n0*exp(B*G*M*m./(y3+R))./exp(B*G*m*M/R);%density at these escape velocities
npn =  n0e*exp(Be*G*M*me./(y3+R))./exp(Be*G*me*M/R);%density at these escape velocities earthlike conditions
npm =  n0e*exp(Be*G*Mm*me./(y3m+Rm))./exp(Be*G*me*Mm/Rm);%density at these escape velocities earthlike conditions moon
npv =  n0e*exp(Be*G*Mc*me./(y3v+Rc))./exp(Be*G*me*Mc/Rc);%density at these escape velocities earthlike conditions ceres
neesc = np.*Pv;%density escaping
neescm = npm.*Pvn;%density escaping moon
neescv = npv.*Pvn;%density escaping ceres
neescn = npn.*Pvn;%density escaping earthlike conditions
cumnp = cumsum(np);
cumnpn = cumsum(npn);
cumnpm = cumsum(npm);
cumnpv = cumsum(npv);
cumne = cumsum(neesc);
cumnen = cumsum(neescn);
cumnem = cumsum(neescm);
cumnev = cumsum(neescv);
perce = cumne(round(sqrt(2*G*M./(R))))/cumnp(round(sqrt(2*G*M./(R))))%percent of all particles at escape velocity
percen = cumnen(round(sqrt(2*G*M./(R))))/cumnpn(round(sqrt(2*G*M./(R))))%percent of all particles at escape velocity earthlike conditions
percem = cumnem(round(sqrt(2*G*Mm./(Rm))))/cumnpm(round(sqrt(2*G*Mm./(Rm))))%percent of all particles at escape velocity earthlike conditions moon
percev = cumnev(round(sqrt(2*G*Mc./(Rc))))/cumnpv(round(sqrt(2*G*Mc./(Rc))))%percent of all particles at escape velocity earthlike conditions ceres