Input, Calling function, Data, Func doesnt work??
显示 更早的评论
So after much re-writing I can't find a solution to why my function doesn't work. I appreciate it's probably written very badly but it results in RA = [empty sym] Dec = [empty sym]. I don't think the variables of the nested functions are being used by the parent function?
Any advice would be greatly appreciated, please and thank you.
function planet_pos()
long = []; inc = []; w = []; a = []; ec = []; m = []; obl = []; ea = [];
v = []; r = []; xec = []; yec = []; zec = []; xeq = []; yeq = []; zeq = [];
year = input('please enter the year');
mon = input('please enter the month as a number');
day = input('please enter the day');
d = 367 * year - (7*(year +((mon + 9)/12)))/4 + (275 * mon)/9 + day - 730530;
planet = input('please enter the planet','s');
function [long,inc,w,a,ec,m,obl]= mercury(d)
long = 48.3313 + 3.24587e-5 * d;
inc = 7.0047 + 5.00e-8 * d;
w = 29.1241 + 1.01444e-5 * d;
a = 0.387098;
ec = 0.205635 + 5.59e-10 * d;
m = 168.6562 + 4.0923344368 * d;
obl = 0.1;
end
function [long,inc,w,a,ec,m,obl]= venus(d)
long = 76.6799 + 2.46590e-5 * d;
inc = 3.3946 + 2.75e-8 * d;
w = 54.8910 + 1.38374e-5 * d;
a = 0.723330;
ec = 0.006773 - 1.302e-9 * d;
m = 48.0052 + 1.6021302244 * d;
obl = 177.4;
end
function [long,inc,w,a,ec,m,obl]= earth(d)
long = 174.8731758 - 0.2410908 * d;
inc = 0 + 0.0130548 * d;
w = 282.9404 + 4.70935e-5;
a = 1.0000010;
ec = 0.00167086 - 1.151e-9 * d;
m = 365.0470 + 0.9856002585 *d;
obl = 23.45;
end
function [long,inc,w,a,ec,m,obl]= mars(d)
long = 49.5574 + 2.11081e-5 * d;
inc = 1.8497 - 1.78E-8 * d;
w = 286.5016 + 2.92961E-5 * d;
a = 1.523688;
ec = 0.093405 + 2.516E-9 * d;
m = 18.6021 + 0.5240207766 * d;
obl = 25.19;
end
function [long,inc,w,a,ec,m,obl]= jupiter (d)
long = 100.4542 + 2.76854e-5 * d;
inc = 1.3030 - 1.557e-7 * d;
w = 273.8777 + 1.64505e-5 * d;
a = 5.20256;
ec = 0.048498 + 4.469e-9 * d;
m = 19.8950 + 0.0830853001 * d;
obl = 3.12;
end
function [long,inc,w,a,ec,m,obl]= saturn(d)
long = 113.6634 + 2.38980e-5 * d;
inc = 2.4886 - 1.081e-7 * d;
w = 339.3939 + 2.97661e-5 * d;
a = 9.55475;
ec = 0.055546 - 9.499e-9 * d;
m = 316.9670 + 0.0334442282 * d;
obl = 26.73;
end
function [long,inc,w,a,ec,m,obl]= uranus(d)
long = 74.0005 + 1.3978e-5 * d;
inc = 0.7733 + 1.9e-8 * d;
w = 96.6612 + 3.0565e-5 * d;
a = 19.18171 - 1.55e-8 * d;
ec = 0.047318 + 7.45e-9 * d;
m = 142.5905 + 0.011725806 * d;
obl = 97.86;
end
function [long,inc,w,a,ec,m,obl]= neptune (d)
long = 131.7806 + 3.0173e-5 * d;
inc = 1.7700 - 2.55e-7 * d;
w = 272.8461 - 6.027e-6 * d;
a = 30.05826 + 3.313e-8 * d;
ec = 0.008606 + 2.15e-9 * d;
m = 260.2471 + 0.005995147 * d;
obl = 29.56;
end
pi20 = vpa(pi);
sind = @(x) sin(x*pi20/180);
asind = @(x) sin(x) * 180/pi20;
cosd = @(x) cos(x*pi20/180);
tand = @(x) tan(x*pi20/180);
atand = @(x) atan(x) * 180/pi20;
atan2d = @(y,x) atan2(y, x) * 180/pi20;
%find eccentric anomaly
ea = m + (180/pi20)* ec * sind(m)* (1 + ec *cosd(m));
keeplooping = true;
while ea > 0.005
E1 = (ea - (180/pi20)* ec *sind(ea)-m)/(1-ec *cosd(ea));
keeplooping;
if ea <= 0.005
keeplooping = false;
end
E1 = ea;
end
%find true anomaly and radial distance
v = 2*atand((sqrt(1+ec))* tand(ea/2));
r = (a *(1-(ec^2))/(1+ec* cosd(v)));
%compute ecliptic rectangular coordinates
x = [];
y = [];
x = cosd(ea) - ec;
y = sind(ea) * sqrt(1 - ec* ec);
r = sqrt(x*x + y*y);
v = atan2d(y,x);
long = v+w;
if long < 0
long = 360 + long;
elseif long > 360
long = long - 360;
end
xec = r * (cosd(long)* cosd(v+w) - sind(long)* sind(v+w)* cosd(inc));
yec = r * (sind(long)* cosd(v+w) - cosd(long)* sind(v+w)* cosd(inc));
zec = r * sind(v+w)* sind(inc);
%compute equitoral coordinates, right ascension and declination then print
xeq = xec;
yeq = yec * cosd(obl) - zec * sind(obl);
zeq = yec * sind(obl) + zec * cosd(obl);
RA = atan2d(yeq, xeq)
eq = sind(obl);
Dec = asind(eq/r)
fprintf('Right Acension is %.2f\n', RA)
fprintf ('Declination is %.2f\n', Dec)
end
9 个评论
You never have code to decode the input planet string and use it to call one of the several planet-specific functions. Just because the variable name planet is assigned in the main function doesn't cause any correlation to occur between that variable and its value and the internal functions of the same string value; the functions must be called by code in the main function.
Without getting into really esoteric things using eval, the ways to do that are to build a CASE statement that selects the function or to make a function handle using the @ operator.
Stephen23
2017-6-3
It could be much simpler to get rid of all of those nested functions and replace them with a simple non-scalar structure, then there is not need for esoteric programming or for switch & case:
planets(1).name = 'Mercury';
planets(1).long = ...
planets(1).inc = ...
...
planets(2).name = 'Venus';
planets(2).long = ...
...
Accessing and vectorizing operations with that data would then be trivial too, for example to get all of the names in one cell array:
names = {planets.name}
dpb
2017-6-3
Maybe it's not enough coffee yet today but I'm not seeing the trick to use that with the variable d as functional, Stephen???
dpb
2017-6-3
...
keeplooping = true;
while ea > 0.005
E1 = (ea - (180/pi20)* ec *sind(ea)-m)/(1-ec *cosd(ea));
keeplooping;
if ea <= 0.005
keeplooping = false;
end
E1 = ea;
end
Above is infinite loop-- ea is computed before the while commences but is never modified inside the loop. In fact, each pass simply resets E1 to the initial value of ea and so nothing effectively ever happens.
Should the condition be on E1 instead, perhaps?
dpb
2017-6-3
Another note on "the way Matlab works" with nested functions--in writing
function [long,inc,w,a,ec,m,obl]= mercury(d)
you've created local variables of all the outputs as well as the input d; these are NOT the variables of the same name in the outer scope but "dummy" variables of the same name with local scope only within the nested function. Thus there's no advantage in writing the functions as nested this way; they could be local functions instead.
To use them as nested functions you would simply write
function mercury
....
end
so all variables are from the parent scope.
With that change, your code would actually function (other than the convergence problem noted above) if you were to use
% evaluate planet function desired
eval(planet)
%find eccentric anomaly
ea = m + (180/pi20)* ec * sind(m)* (1 + ec *cosd(m));
...
Normally eval is frowned upon, but this invocation is simple enough as to not be difficult to debug the command and does do what you're looking for in how you've written the overall function.
I'd probably restructure it differently, but this isn't too bad, imo.
Would probably use local functions and array of function handles or the like, but can see the above given there is no complex command being built for eval.
dpb
2017-6-3
...
pi20 = vpa(pi);
sind = @(x) sin(x*pi20/180);
asind = @(x) sin(x) * 180/pi20;
cosd = @(x) cos(x*pi20/180);
tand = @(x) tan(x*pi20/180);
atand = @(x) atan(x) * 180/pi20;
atan2d = @(y,x) atan2(y, x) * 180/pi20;
Why all the above, Matlab has builtin trig functions of the same names for inputs in degrees, why not use them?
As well as don't need symbolic to return pi to machine precision; it's also builtin function.
Perhaps the above came from Octave or somewhere else that might not have had the degree counterparts, or just unaware of them?
BTW, in the above asind appears to be missing an a on the RH side if you do use these...
I am not convinced that nested or even local functions are required. Why not simply use a structure:
S.mercury.long = ...
S.mercury.inc = ...
S.mercury....
...
S.venus.long = ...
S.venus.inc = ...
... etc
planet = input('please enter the planet','s');
data = S.(planet);
and all of the values are there:
data.long
data.inc
data... etc
It would even be possible to push these into separate varaibles (although I woud not bother):
long = data.long
inc = data.inc
etc
Why make this so complicated?
dpb
2017-6-3
I'm still not seeing that the structure solves the computation question, Stephen?
Granted, it has many more level than needed; there could be simply one function that populates the result with a set--oh, I guess that's what you're saying...gotcha'.
Cat
2017-6-6
回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Mathematics 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
