EarthToMars.zip_EarthtoMars_lambert matlab_orbit optimization_地

% LAMBERT Lambert-targeter for ballistic flights % (Izzo, and Lancaster, Blanchard & Gooding) % % Usage: % [V1, V2, extremal_distances, exitflag] = lambert(r1, r2, tf, m, GM_central) % % Dimensions: % r1, r2 -> [1x3] % V1, V2 -> [1x3] % extremal_distances -> [1x2] % tf, m -> [1x1] % GM_central -> [1x1] % % This function solves any Lambert problem *robustly*. It uses two separate % solvers; the first one tried is a new and unpublished algorithm developed % by Dr. D. Izzo from the European Space Agency [1]. This version is extremely % fast, but especially for larger [m] it still fails quite frequently. In such % cases, a MUCH more robust algorithm is started (the one by Lancaster & % Blancard [2], with modifcations, initial values and other improvements by % R.Gooding [3]), which is a lot slower partly because of its robustness. % % INPUT ARGUMENTS: % ====================================================================== % name units description % ====================================================================== % r1, r1 [km] position vectors of the two terminal points. % tf [days] time of flight to solve for % m [-] specifies the number of complete orbits to complete % (should be an integer) % GM_central [km3/s2] std. grav. parameter (G愚 = mu) of the central body % % OUTPUT ARGUMENTS: % ====================================================================== % name units description % ====================================================================== % V1, V2 [km/s] terminal velocities at the end-points % extremal_distances [km] minimum(1) and maximum(2) distance of the % spacecraft to the central body. % exitflag [-] Integer containing information on why the % routine terminated. A value of +1 indicates % success; a normal exit. A value of -1 % indicates that the given problem has no % solution and cannot be solved. A value of -2 % indicates that both algorithms failed to find % a solution. This should never occur since % these problems are well-defined, and at the % very least it can be determined that the % problem has no solution. Nevertheless, it % still occurs sometimes for accidental % erroneous input, so it provides a basic % mechanism to check any application using this % algorithm. % % This routine can be compiled to increase its speed by a factor of about % 50, which is certainly advisable when the complete application requires % a great number of Lambert problems to be solved. The entire routine is % written in embedded MATLAB, so it can be compiled with the emlmex() % function. To do this, make sure MATLAB's current directory is equal to % where this file is located. Then, copy-paste and execute the following % commands to the command window: % % example_input = {... % [0.0, 0.0, 0.0], ...% r1vec % [0.0, 0.0, 0.0], ...% r2vec % 0.0, ... % tf % 0.0, ... % m % 0.0}; % muC % emlmex -eg example_input lambert.m % % This is of course assuming your compiler is configured correctly. See the % docs on emlmex() on how to do that. % % References: % % [1] Izzo, D. ESA Advanced Concepts team. Code used available in MGA.M, on % http://www.esa.int/gsp/ACT/inf/op/globopt.htm. Last retreived Nov, 2009. % [2] Lancaster, E.R. and Blanchard, R.C. "A unified form of Lambert's theorem." % NASA technical note TN D-5368,1969. % [3] Gooding, R.H. "A procedure for the solution of Lambert's orbital boundary-value % problem. Celestial Mechanics and Dynamical Astronomy, 48:145�165,1990. % % See also lambert_low_ExpoSins. % Author: % Name : Rody P.S. Oldenhuis % E-mail : oldenhuis@dds.nl / oldenhuis@gmail.com % Affiliation: Delft University of Technology % % please report any bugs or suggestions to oldnhuis@dds.nl. % 毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽� % Izzo's version: % Very fast, but not very robust for more complicated cases % 毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽� function [V1, V2, extremal_distances, exitflag] = lambert(r1vec, r2vec, tf, m, muC)%#eml % original documentation: % 毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽毽� % % This routine implements a new algorithm that solves Lambert's problem. The % algorithm has two major characteristics that makes it favorable to other % existing ones. % % 1) It describes the generic orbit solution of the boundary condition % problem through the variable X=log(1+cos(alpha/2)). By doing so the % graph of the time of flight become defined in the entire real axis and % resembles a straight line. Convergence is granted within few iterations % for all the possible geometries (except, of course, when the transfer % angle is zero). When multiple revolutions are considered the variable is % X=tan(cos(alpha/2)*pi/2). % % 2) Once the orbit has been determined in the plane, this routine % evaluates the velocity vectors at the two points in a way that is not % singular for the transfer angle approaching to pi (Lagrange coefficient % based methods are numerically not well suited for this purpose). % % As a result Lambert's problem is solved (with multiple revolutions % being accounted for) with the same computational effort for all % possible geometries. The case of near 180 transfers is also solved % efficiently. % % We note here that even when the transfer angle is exactly equal to pi % the algorithm does solve the problem in the plane (it finds X), but it % is not able to evaluate the plane in which the orbit lies. A solution % to this would be to provide the direction of the plane containing the % transfer orbit from outside. This has not been implemented in this % routine since such a direction would depend on which application the % transfer is going to be used in. % % please report bugs to dario.izzo@esa.int % % adjusted documentation: % 毽毽毽毽毽毽毽毽毽毽毽� % % By default, the short-way solution is computed. The long way solution % may be requested by giving a negative value to the corresponding % time-of-flight [tf]. % % For problems with |m| > 0, there are generally two solutions. By % default, the right branch solution will be returned. The left branch % may be requested by giving a negative value to the corresponding % number of complete revolutions [m]. % Authors % .戡�.戡�.戡�.戡�.戡�.戡�.戡�.戡�.戡�.戡�.戡�.戡�.戡�.戡�.戡�. % Name : Dr. Dario Izzo % E-mail : dario.izzo@esa.int % Affiliation: ESA / Advanced Concepts Team (ACT) % Made readible and optimized for speed by Rody P.S. Oldenhuis % Code available in MGA.M on http://www.esa.int/gsp/ACT/inf/op/globopt.htm % last edited 12/Dec/2009 % ADJUSTED FOR EML-COMPILATION 24/Dec/2009 % initial values tol = 1e-12; bad = false; days = 86400; % work with non-dimensional units r1 = sqrt(r1vec*r1vec.'); r1vec = r1vec/r1; V = sqrt(muC/r1); r2vec = r2vec/r1; T = r1/V; tf = tf*days/T; % also transform to seconds % relevant geometry parameters (non dimensional) mr2vec = sqrt(r2vec*r2vec.'); % make 100% sure it's in (-1 <= dth <= +1) dth = acos( max(-1, min(1, (r1vec*r2vec.')/mr2vec)) ); % dec