Code covered by the BSD License  

Highlights from
Gravity

image thumbnail
from Gravity by link2358 noname
Computes position of several bodies under the action of gravitatory forces.

secant1(inibod,endbod,inid,init,dt2,tol,tol2,miter,inig)
function [alpha1, iniv1, error, kiter] = secant1(inibod,endbod,inid,init,dt2,tol,tol2,miter,inig)

%once gravity 2 has been run, and given the entry parameters (starting body, time ,
%target body) , it computes the necessary angle to impact the target
%body with a starting velocity equal to the escape velocity
%IT MAY NOT HAVE A SOLUTION FOR CERTAIN INITIAL CONDITIONS

%inibod: starting body
%endbod: target body
%inid: distance from the initial body from where the rocket is launched
%dt2: time step size
%init: launching time
%inidir: launching direction [cos(alpha; sin(alpha)]
%tol and tol2: tolerances for numerical iterative processes
%miter: maximum number of iterations iterative processes

close all
dbstop if error
global g mass ndime  
   
if ndime == 3
    disp('Not ready for 3 dimensional case');
    stop
end


alpha0 = inig;
alpha1 = alpha0-0.01;

x0 = alpha0;
x1 = alpha1;

ve = sqrt(2*g*mass(inibod)/inid);

[mdist] = mindist(inibod,inid,dt2,init,[cos(x0);sin(x0)],ve,tol2,miter,Inf);
fx0 = mdist(endbod);
[mdist] = mindist(inibod,inid,dt2,init,[cos(x1);sin(x1)],ve,tol2,miter,Inf);
fx1 = mdist(endbod);


error = 1;
kiter = 0;

while error > tol && kiter < miter
    kiter = kiter +1;
    grad1 = (fx1-fx0)/(x1-x0);

    %grad1 = grad_bo([fx1,fx0,fx_1],[x1';x0';x_1']);
    xk = x1 - fx1/grad1; %avanza en la direccin del gradiente
    [mdist] = mindist(inibod,inid,dt2,init,[cos(xk);sin(xk)],ve,tol2,miter,Inf);
    fxk = mdist(endbod);
    
    error = abs(fxk);
    x0 = x1;
    fx0 = fx1;
    x1 = xk;
    fx1 = fxk;
end

if x1 > 2*pi
    x1 = rem(x1,2*pi);
elseif x1 < 0
    x1 = rem(x1,2*pi)+2*pi;
end
alpha1 = x1;
iniv1 = ve;

[mdist] = mindist(inibod,inid,dt2,init,[cos(x1);sin(x1)],ve,tol2,miter,1);

Contact us at files@mathworks.com