@Carl Love 

I mean an argument for the definition where the label should be positioned in the 3d plot. Previous posts suggest something like ...,[1,3],... in a 2d plot. I tried something like this in the above 3d plot but it does not work.

@Carl Love 

> with(linalg);

> with(plots);

> with(plottools);

> P1 := matrix([[cos(omega[j]), -sin(omega[j]), 0], [sin(omega[j]), cos(omega[j]), 0], [0, 0, 1]]); P2 := matrix([[1, 0, 0], [0, cos(i[j]), -sin(i[j])], [0, sin(i[j]), cos(i[j])]]); P3 := matrix([[cos(Omega[j]), -sin(Omega[j]), 0], [sin(Omega[j]), cos(Omega[j]), 0], [0, 0, 1]]);

> A:=matrix([[a[j]*(cos(E[j])-e[j])],[a[j]*sqrt(1-e[j]^2)*sin(E[j])],[0]]);

> R:=multiply(P3,P2,P1);

> B:=multiply(R,A);

> a := [.38709893, .72333199, 1.00000011, 1.52366231, 5.20336301, 9.53707032, 19.19126393, 30.06896348, 39.48168677, 1/0.5454e-2, 268.2509283, 532.7838156, 300];

> e := [.20563069, 0.677323e-2, 0.1671022e-1, 0.9341233e-1, 0.4839266e-1, 0.5415060e-1, 0.4716771e-1, 0.858587e-2, .24880766, .994920, .7005635, .8570973, .1];

> i := [7.00487, 3.39471, 0.5e-4, 1.85061, 1.30530, 2.48446, .76986, 1.76917, 17.14175, 89.5328, 24.01830, 11.92859, 10];

> Omega := [48.33167, 76.68069, -11.26064, 49.57854, 100.55615, 113.71504, 74.22988, 131.72169, 110.30347, 282.1476, 90.88303, 144.53190, 45];

> omega := [77.45645, 131.53298, 102.94719, 336.04084, 14.75385, 92.43194, 170.96424, 44.97135, 224.06676, 130.8038, 293.03200, 311.18311, 150];

> i := map(x→ convert(x, units, deg, rad) end proc, i);

> Omega := map(x→ convert(x, units, deg, rad) end proc, Omega);

> omega := map(x→ convert(x, units, deg, rad) end proc, omega);

> for j to 13 do omega[j] := arcsin(sin(omega[j]-Omega[j])/sin(arccos(sin(i[j])*cos(omega[j]-Omega[j])))) end do;

> x := array(1 .. 13);

> y := array(1 .. 13);

> z := array(1 .. 13);

> for j to 13 do x[j] := B[1, 1]; y[j] := B[2, 1]; z[j] := B[3, 1] end do;

> Sun := array([1]);

> Inner := array(1 .. 18); for j to 9 do Colors := [black, green, blue, red, black, yellow, green, violet, brown, aquamarine, black, black, red]; Linestyle := [solid, solid, solid, solid, solid, solid, solid, solid, solid, longdash, solid, solid, longdash]; Inner[j] := spacecurve([subs(E[j] = E, x[j]), subs(E[j] = E, y[j]), subs(E[j] = E, z[j])], E = 0 .. 2*Pi, color = Colors[j], linestyle = Linestyle[j]) end do;

> Comet := array([1, 2]); if j = 10 then Colors := [aquamarine]; Linestyle := [longdash]; Comet[1] := spacecurve([subs(E[j] = E, x[j]), subs(E[j] = E, y[j]), subs(E[j] = E, z[j])], E = 0 .. 2*Pi, color = Colors[j], linestyle = Linestyle[j]) end if;

> Oort := array(1 .. 6); for j from 11 to 13 do Colors := [black, black, red]; Linestyle := [solid, solid, longdash]; Inner[j] := spacecurve([subs(E[j] = E, x[j]), subs(E[j] = E, y[j]), subs(E[j] = E, z[j])], E = 0 .. 2*Pi, color = Colors[j], linestyle = Linestyle[j]) end do;

> Sun[1] := plot3d(0.1e-1, 0 .. 2*Pi, 0 .. Pi, style = PATCHNOGRID, coords = spherical, color = red);

> Inner[10] := textplot3d([subs(E[1] = 0, x[1]), subs(E[1] = 0, y[1]), subs(E[1] = 0, z[1]), "Mercury"]); Inner[11] := textplot3d([subs(E[2] = 0, x[2]), subs(E[2] = 0, y[2]), subs(E[2] = 0, z[2]), "Venus"]); Inner[12] := textplot3d([subs(E[3] = 0, x[3]), subs(E[3] = 0, y[3]), subs(E[3] = 0, z[3]), "Earth"]); Inner[13] := textplot3d([subs(E[4] = 0, x[4]), subs(E[4] = 0, y[4]), subs(E[4] = 0, z[4]), "Mars"]); Inner[14] := textplot3d([subs(E[5] = 0, x[5]), subs(E[5] = 0, y[5]), subs(E[5] = 0, z[5]), "Jupiter"]); Inner[15] := textplot3d([subs(E[6] = 0, x[6]), subs(E[6] = 0, y[6]), subs(E[6] = 0, z[6]), "Saturn"]); Inner[16] := textplot3d([subs(E[7] = 0, x[7]), subs(E[7] = 0, y[7]), subs(E[7] = 0, z[7]), "Uranus"]); Inner[17] := textplot3d([subs(E[8] = 0, x[8]), subs(E[8] = 0, y[8]), subs(E[8] = 0, z[8]), "Neptune"]); Inner[18] := textplot3d([subs(E[9] = 0, x[9]), subs(E[9] = 0, y[9]), subs(E[9] = 0, z[9]), "Pluto"]); Comet[2] := textplot3d([subs(E[10] = 0, x[10]), subs(E[10] = 0, y[10]), subs(E[10] = 0, z[10]), Hale-Bopp]); Oort[4] := textplot3d([subs(E[11] = 0, x[11]), subs(E[11] = 0, y[11]), subs(E[11] = 0, z[11]), "2012 VP113"]); Oort[5] := textplot3d([subs(E[12] = 0, x[12]), subs(E[12] = 0, y[12]), subs(E[12] = 0, z[12]), "Sedna"]); Oort[6] := textplot3d([subs(E[13] = 0, x[13]), subs(E[13] = 0, y[13]), subs(E[13] = 0, z[13]), "Planet X ?", color = red]);

> Sun1 := convert(Sun, 'list');
> Inner1 := convert(Inner, 'list');
> Comet1 := convert(Comet, 'list');
> Oort1 := convert(Oort, 'list');
> display(Sun1, Inner1, Comet1, Oort1, scaling = CONSTRAINED);

 

The first error message appears after the if-condition. Can you tell show me where I am making a mistake? Beware: when copy paste the code from Maple to Word and from Word in here the colons at the end of lines have changed into semi-colons. Hope this is no problem in executing the code despite the lines being in the same ">..." .

@Axel Vogt 

me too

@lt 

Moreover, trying to color fill the asteroid belt:

     > MainBelt := seq([pic[15] .. pic[16]], color = brown);
     > display({pic1, MainBelt}, scaling = CONSTRAINED);
seems not to work

I added two more circles for the borders of the asteroid belt, changed the numbers from 9 to 11, added two more orbital parameters(just 2.3 and 3.2 AU with all other values zero). But now the circles for Neptune and Pluto are missing although their label is shown.

Orbital parameters (Semi-major axis [AU], Eccentricity, Inclination[deg], Longitude of ascending node [deg], Longitude of Perihelion[deg]):
    > a := [.38709893, .72333199, 1.00000011, 1.52366231, 2.2, 3.2, 5.20336301, 9.53707032, 19.19126393, 30.06896348, 39.48168677];

    > e := [.20563069, 0.677323e-2, 0.1671022e-1, 0.9341233e-1, 0, 0, 0.4839266e-1, 0.5415060e-1, 0.4716771e-1, 0.858587e-2, .24880766];

    > i := [7.00487, 3.39471, 0.5e-4, 1.85061, 0, 0, 1.30530, 2.48446, .76986, 1.76917, 17.14175];

    > Omega := [48.33167, 76.68069, -11.26064, 49.57854, 0, 0, 100.55615, 113.71504, 74.22988, 131.72169, 110.30347];

    > omega := [77.45645, 131.53298, 102.94719, 336.04084, 0, 0, 14.75385, 92.43194, 170.96424, 44.97135, 224.06676];

Degree to Radian conversion:
    > i := map(proc (x) options operator, arrow; convert(x, units, deg, rad) end proc, i);

    > Omega := map(proc (x) options operator, arrow; convert(x, units, deg, rad) end proc, Omega);

    > omega := map(proc (x) options operator, arrow; convert(x, units, deg, rad) end proc, omega);

    > for j to 11 do

    omega[j] := arcsin(sin(omega[j]-Omega[j])/sin(arccos(sin(i[j])*cos(omega[j]-Omega[j])))) end do;

Parameter equations:
    > x := array(1 .. 11); y := array(1 .. 11); z := array(1 .. 11);

    > for j to 11 do

    x[j] := B[1, 1]; y[j] := B[2, 1]; z[j] := B[3, 1] end do;

    > pic := array(1 .. 21);

    > for j to 11 do

    Colors := [black, green, blue, red, brown, brown, black, gold, violet, green, brown];

    >pic[j] := spacecurve([subs(E[j] = E, x[j]), subs(E[j] = E, y[j]), subs(E[j] = E, z[j])], E = 0 .. 2*Pi, color = Colors[j]) end do;

    > pic[10] := plot3d(0.1e-1, 0 .. 2*Pi, 0 .. Pi, style = PATCHNOGRID, coords = spherical, color = red);

    > pic[11] := textplot3d([subs(E[1] = 0, x[1]), subs(E[1] = 0, y[1]), subs(E[1] = 0, z[1]), "Mercury"]);

    > pic[12] := textplot3d([subs(E[2] = 0, x[2]), subs(E[2] = 0, y[2]), subs(E[2] = 0, z[2]), "Venus"]);

    > pic[13] := textplot3d([subs(E[3] = 0, x[3]), subs(E[3] = 0, y[3]), subs(E[3] = 0, z[3]), "Earth"]);

    > pic[14] := textplot3d([subs(E[4] = 0, x[4]), subs(E[4] = 0, y[4]), subs(E[4] = 0, z[4]), "Mars"]);

    > pic[15] := textplot3d([subs(E[5] = 0, x[5]), subs(E[5] = 0, y[5]), subs(E[5] = 0, z[5]), "MABin"]);

    > pic[16] := textplot3d([subs(E[6] = 0, x[6]), subs(E[6] = 0, y[6]), subs(E[6] = 0, z[6]), "MABout"]);

    > pic[17] := textplot3d([subs(E[7] = 0, x[7]), subs(E[7] = 0, y[7]), subs(E[7] = 0, z[7]), "Jupiter"]);

    > pic[18] := textplot3d([subs(E[8] = 0, x[8]), subs(E[8] = 0, y[8]), subs(E[8] = 0, z[8]), "Saturn"]);

    > pic[19] := textplot3d([subs(E[9] = 0, x[9]), subs(E[9] = 0, y[9]), subs(E[9] = 0, z[9]), "Uranus"]);

    > pic[20] := textplot3d([subs(E[10] = 0, x[10]), subs(E[10] = 0, y[10]), subs(E[10] = 0, z[10]), "Neptune"]);

    > pic[21] := textplot3d([subs(E[11] = 0, x[11]), subs(E[11] = 0, y[11]), subs(E[11] = 0, z[11]), "Pluto"]);

    > pic1 := convert(pic, 'list');

    > display(pic1, scaling = CONSTRAINED);

@tomleslie 

Yes exactly. I need one single, warped suface going through all the circles! Not plane surfaces filling each circle. Sorry for not specifying this enough in the entry. In the header I also could not write the whole problem and the current one with the circles seemed a good entry to me. Sorry for the confusion!

Concerning the professional post doc: You are telling me your opinion - btw, to despise him is not very polite. He was telling me his experience. Many others have also told me that physicists in academia are not known for their programming skills. Informaticians and software engineers are usually better. Notice the "usually", of course - which is enough of differentiated thinking in this context. That's all I am writing and the others were telling - not more, not less. Let's leave this topic and get back to the warped surface.

The problem must lie at the point where the calculated inclination, radius, and phase values are put into the function of the surface points.

P.S. This "Project Euler" seems impressive! For sure this is extremely useful for refreshing/exercising programming skills.

Thank you very much for your effort to help me – all of you. Just a short notice before continuing. Phrases like your reply-title above are unpolite and discouraging. First, physicists are not usually known for their ability to make beautiful and wonderful codes compared to informaticians and software engineers – especially students and rookies as I am in Maple. Second, from a professional post-doc in theoretical physics and himself having a professional education and work experience in industrial software design (therefore knowing how a good code is being built up) I know that in academia codes are usually not produced for beauty-reasons and have to work just for this one task and for the next paper you can throw them into the dustbin. So, please, not such phrases as above anymore.

I post here a picture of what I want to achieve. You see an R=30 and 50 inner and outer circle with a blue and yellow circle in-between. These are 4 inclination points when looking in one direction from the center. In orange you see the 100 spline-interpolated inclination points. The more points there are the smother the final surface will be. Now make this interpolation for all angles from 0 to 359° and plot a surface through all inclination points. The beginning of such a surface is indicated on the right side. In the end, I will see if I want to have circles indicating specific radii and axes or color coding, e.g. from blue (minus inclinations) to yellow (plus inclinations) for clarification.warped_plane.docx

What is wrong with Mapleprimes? I can't upload Wordfiles, jpg, etc. The connection always breaks down. This is not helpfull at all!

@tomleslie 
now it shows me only semi-circles. I'm sure this can be solved with some kind of range-altering. But what I tried to say with my interpolation approach was that I want to have one single surface that goes through all circles, not a plane surface filling each circle. I need one single surface that visualises the inclination depending on orbital radius and angle seen from the center of the circles. What I have achieved so far is to interpolate between arbitrary inclination and phases as a function of radius. The only thing left to do is to put these values into your equations and plot them with plot3d.
     > R := [30, 32.5, 37.5, 42.5, 47.5, 50];
     > incl := [0, .5, 1, -2, 2*evalf(Pi)-1.75, 0];
     > pointplot(R, incl, color = [grey, red, blue, green, black, grey], labels = ["radius", "incl"]);
     > NewR := [seq(30+.2*i, i = 0 .. 100)];
     > Newincl := ArrayInterpolation(R, incl, NewR, method = spline);
     > pointplot(NewR, Newincl, labels = ["radius", "incl"]);
     > phases := [0, (1/12)*evalf(Pi), 3*evalf(Pi)*(1/4), (1/4)*evalf(Pi)+.2, evalf(Pi)/(2.5), 0];
     > pointplot(R, phases, labels = ["radius", "phase"]);
     > Newphases := ArrayInterpolation(R, phases, NewR, method = spline);
     > pointplot(NewR, Newphases, labels = ["radius", "phase"]);
     > f(i,t):=[NewR[i]*cos(t),NewR[i]*sin(t),Newincl[i]*cos(t-Newphases[i])];
     > plot3d(f, t=0..2*Pi, view=[-50..50,.50..50,-15--15]);

Everything works except for the last plot3d. Although f depends on i I did not specify this because the values were already calculated in the inerpolations using i=0..100, so should not be specified anymore in the last plot.

@Carl Love 

> R := [30, 32.5, 37.5, 42.5, 47.5, 50];
> incl := [0, .5, 1, -2, 2*Pi-1.75, 0];
> pointplot(R, incl, color = [grey, red, blue, green, black, grey], labels = ["Radius", "incl"]);

> NewR := [seq(30+.2*i, i = 0 .. 100)];
> Newincl := ArrayInterpolation(R, incl, NewR, method = spline);
Error, (in CurveFitting:-ArrayInterpolation) invalid input: dependent coordinates must be of type extended_numeric
> pointplot(NewR, Newincl, labels = ["Radius", "incl"]);

@tomleslie

Absolutely. One meta-problem is that my posts seem not to be accepted all the time and I am thrown out and have to loggin again. Not being able to copy my code here as in any other program/site is the second meta-problem. The third is that I did not see that my initial post could be misunderstood. Thus, my clarification of a surface between the circles below:
This is my new approach by interpolation:
I define the beginning and end of my belt at 30 AU and 50 AU and the inclination of the objects at each distance. A pointplot visualises this.

     > R := [30, 32.5, 37.5, 42.5, 47.5, 50];
     > incl := [0, .5, 1, -2, 5, 0];
     > pointplot(R, incl, color = [grey, red, blue, green, black, grey]);

Afterwards, I use the interpolation approach to calculate the inclinations for increases in radius of 0.2 AU. Nice pointplot again.

     > NewR := [seq(30+.2*i, i = 0 .. 100)];
     > Newincl := ArrayInterpolation(R, incl, NewR, method = spline);
     > pointplot(NewR, Newincl);

As I now have enough radius and inclination values to create a smoth surface over my Solar system I used Tom's equations to plot the inclined circles (by the way, for simplification, the orbits of these objects can be assumed as circular as a first approximation; but there are more eccentric object families out there of course). The phases of the given objects are also specified, as well as the interpolated new phases inbetween them. Putting all this into Tom's equations should give me a nice surface plot with varying inclination of radius R and angle t.         

     > t = 0 .. 2*Pi;
     > phases := [0, (1/12)*Pi, (1/4)*Pi, (1/8)*Pi, (1/2)*Pi, 0];
     > Newphases := [seq(0+2*Pi*(1/100), i = 0 .. 100)];
     > f(i,t):=[NewR[i]*cos(t),NewR[i]*sin(t), Newincl[i]*cos(t-Newphases[i])];
     > plot3d([NewR[i]*cos(t), NewR[i]*sin(t), Newincl[i]*cos(t-Newphases[i])], -50 .. 50, -50 .. 50);
     Warning, unable to evaluate the function to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct

What am I missing?

@Carl Love 

While Kitonum4065's code only showed me the simulation at t=0 but nothing more, Tom's code gave me this error message:

> with(plots); ampl := [2.0, 1.1, 1.2, 1.4]; incl := [0., .1, .2, .3]; phases := [0, (1/8)*Pi, (1/4)*Pi, (1/2)*Pi]; colors := [black, blue, red, green]; ps := [seq*(implicitplot3d*(z = incl[j]*r*cos(theta-phases[j]), r = 0 .. ampl[j], theta = -Pi .. Pi, z = -1 .. 1, coords = cylindrical, color = colors[j], style = surface, numpoints = 10000, transparency = .5), j = 1 .. 4)]; display(ps);
Error, (in plots:-display) expecting plot structures but received: [seq*(plots:-implicitplot3d*(z = [0., .1, .2, .3][j]*r*cos(theta-[0, (1/8)*Pi, (1/4)*Pi, (1/2)*Pi][j]), r = 0 .. [2.0, 1.1, 1.2, 1.4][j], theta = -Pi .. Pi, z = -1 .. 1, coords = cylindrical, color = [black, blue, red, green][j], style = surface, numpoints = 10000, transparency = .5), j = 1 .. 4)]

@tomleslie

It should somehow work like this: inside the seq() command the x,y, and z-coordinates of our circles are calculated which are plotted with the overal display-command. I suggest a surface-plotting command inside the display-brackets that use the radius information and z-coordinate of each circle in our 3D plot to interpolate the z-values between the circles. Those spline-interpolations might be useful.

Excuse my long absence...exams.
Thank you all very much for your suggestions. Unfortunately, Tom, your code only gives me en error message instead of a graphic or something. I am using Maple13, I guess you have a more recent version? Based on your code
    with(plots):
    ampl:=[1.0, 1.1, 1.2, 1.4]:
    incl:=[0.0, 0.1, 0.2, 0.3]:
    phases:=[0, Pi/8, Pi/4, Pi/2]:
    colors:=[black, blue, red, green]:
    display(  seq( spacecurve( [ ampl[j]*cos(t),
                                             ampl[j]*sin(t),
                                             incl[j]*cos(t-phases[j])
                                          ],
                                          t=-Pi..Pi,
                                          color=colors[j],
                                          view=[-1.5..1.5,-1.5..1.5,-1..1]
                                       ),
                         j=1..4
                       )
               );
I am trying to plot a warped plane going through these differently inclined circles. Unfortunately, I do not succeed. Imagine we are looking down onto our cicrles and we cut them in half. Now we look from the side and in a 2D x-y-plot we see where our circles are located - careful: the 2D x-y-axes shall not have the same orientation as in the 3D x-y-z-illustration! Ok, we now have our 2D dot-locations where our circles are located depending on the direction we are looking at from the center of our circles. In one direction a circle may have an y-coordinate of 2.34 and when I turn my had 180° it has the y-coordinate -2.34 . My idea now was to use a plot command that interpolates a curve through all my f(x)=y-coordinates begining from and ending on arbitrary x values. E.g. our circles are located at x-values 1.45, 2.3, 4,27 and on y-values -0.5, 1.5, 0.3. I want my interpolation to start and end at x=1.5 .. 4.5 . Depending on where we look from the center outwards the interpolation has to go through the current y-values...or in this 3D case our z-values. So I want the interpolation command to go through 0 to 360° and afterwards plot a surface through the circles as well as the interpolated values. I guess this could have something to do with CurveFitting[ArrayInterpolation]?

Thank you very much! This is very useful! Where would I have to include plot3d(...) in order to show a surface going through all these circles?