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

[dpb]

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]

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]

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]

...
  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]

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]

...
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...

[Stephen23]

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

And then accessing the required data is simple using dynamic fieldnames:

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]

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]

thank you so much, definitely gives me some direction now.. very much appreciated!!