Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- Vinf=[2.7 4 6 8];
- mu_m=42828;
- R_msoi=.576e6;
- [afas,efas,pfas]=parasint(Vinf);
- an_vera0=zeros(1,length(Vinf));
- Vr0=zeros(1,length(Vinf));
- Vt0=zeros(1,length(Vinf));
- for j=1:length(Vinf)
- c_tstar0=(1/efas(j))*((pfas(j)/R_msoi)-1);
- an_vera0(j)=-acos(c_tstar0);
- s_tstar0=sin(an_vera0(j));
- Vr0(j)=sqrt(mu_m/pfas(j))*efas(j)*s_tstar0;
- Vt0(j)=sqrt(mu_m/pfas(j))*(1+efas(j)*c_tstar0);
- end
- ti=time(an_vera0,efas,afas,mu_m);
- for n=1:length(ti)
- I=(0:3600:ti(n));
- options=odeset('Vectorized','on','RelTol',1e-8,'AbsTol',1e-9);
- [t,Y]=ode45(@(t,y) eqMotoCur(t,y),I,y0,options);
- end
- function [afas,efas,pfas]=parasint(Vinf)
- mu_m=42828;
- R_m=3396.2;
- Rp0=R_m+400;
- Rpf=R_m+50;
- a0=zeros(1,length(Vinf));
- e0=zeros(1,length(Vinf));
- p0=zeros(1,length(Vinf));
- afas=zeros(1,length(Vinf));
- efas=zeros(1,length(Vinf));
- pfas=zeros(1,length(Vinf));
- for j=1:length(Vinf)
- a0(j)=(-mu_m)/(Vinf(j)^2);
- e0(j)=1-Rp0/a0(j);
- p0(j)=a0(j)*(1-e0(j)^2);
- c(1)=p0(j);
- c(2)=Rpf*(1-e0(j)^2);
- c(3)=Rpf*(1-e0(j)^2)-p0(j);
- Efas=roots(c);
- ind=(Efas>1 & isreal(Efas));
- efas(j)=Efas(ind);
- afas(j)=Rpf/(1-efas(j));
- pfas(j)=afas(j)*(1-efas(j)^2);
- end
- function [ti] = time(an_vera0,efas,afas,mu_m)
- ti=zeros(1,length(an_vera0));
- for j=1:length(an_vera0)
- F=2*atanh(sqrt((efas(j)-1)/(efas(j)+1))*tan(an_vera0(j)/2));
- T=sqrt((-afas(j))^3/mu_m)*(efas(j)*sinh(F) - F);
- ti(j)=-2*T;
- end
- function [dydt] = eqMotoCur(t,y)
- R_m=3396.2;
- mu_m=42828;
- Cd=2.2;
- S=7.065e-6;
- m=3699.046;
- wE=0.24117/R_m;
- if abs(y(1))>R_m+250 && y(3)>0
- dydt(1:4)=0;
- else
- dydt=zeros(4,1);
- dydt(1)=y(3);
- dydt(2)=y(4)/y(1);
- dydt(3)=(-mu_m/(y(1)^2))+((y(4)^2)/y(1))-(.5*Cd)*(S/m)*...
- (dens(y(1)-R_m))*(sqrt((y(3)^2)+(y(4)-(wE*y(1)))^2))*y(3);
- dydt(4)=-(y(4)*(y(3)/y(1)))-(.5*Cd)*(S/m)*...
- (dens(y(1)-R_m))*(sqrt((y(3)^2)+(y(4)-(wE*y(1)))^2))*...
- (y(4)-(wE*y(1)));
- end
- end
- function [rho] = dens(z)
- load Density.mat h d
- d=d*10^9;
- if any(h)==abs(z)
- j=h==abs(z);
- rho=d(j);
- elseif abs(z)<250
- c_tot=(find(h>=abs(z)));
- c=c_tot(1);
- H=-(h(c)-h(c-1))/(log(d(c)/d(c-1)));
- rho=d(c-1)*exp(-(abs(z)-h(c-1))/H);
- else
- rho=0;
- end
- end
- Unable to perform assignment because the size of the left side is
- 4-by-2 and the size of the right side is 4-by-5.
- Error in ode45 (line 488)
- yout(:,idx) = yout_new;
- Error in MainCur (line 59)
- [t,Y]=ode45(@(t,y) eqMotoCur(t,y),I,y0,options);
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement