1
$\begingroup$

Thanks to Henrik Schumacher for advice to simplify the sample code:

I would be grateful for any pointers on the following issue.

A line from an external point intersects the surface at 2 points (near side and far side) The line can be defined as emanating from the origin {0,0,0} with a horizontal angle L and a vertical angle B (or in Cartesian form) I have tried several approaches to identify the distance along the line for each intersection with the 3D surface.

Obviously given any point on the surface this is straight forward and one might think that it would also be straightforward to look for two points on the same line (i.e. same L and B) and thus calculate the distance along the line to each intersection , but this has proved intractable as there are rarely exact matches and more data points creates confusion.

I have tried several approaches including meshing but to no avail, they just seem to give another set of points that do not seem to be triangles that share edges, but that may be my misunderstanding of what a mesh does.

I note that If there is a single equation rather than a parametric set of them then NSolve or Solve can be used.

The example below includes enough code to generate a sample surface the parameters are u and v and the x ,y and z values are given by rather involved functions of which these are simplified examples.


    ClearAll;

    fx[i_, j_, v_, h_] := 8.5 - j (1 - Sin[v]) Cos[v] Cos[h];

    fy[i_, j_, v_, h_] := i (1 - .2*Sin[v]) Cos[v] Sin[h];

    fz[k_, v_] := k (1 - Sin[v]);




    Manipulate[
     FB = ParametricPlot3D[{
        fy[i, j, v, h],
        fx[i, j, v, h],
        fz[k, v]
        },
       {v, 0, 2 Pi},
       {h, -Pi/2, Pi/2},

       Mesh -> Automatic
       ],

     {{i, 35, "y axis (i)"}, 5, 70, Appearance -> "Labeled"},
     {{j, 35, "x axis (j)"}, 5, 70, Appearance -> "Labeled"},
     {{k, 50, "z axis (k)"}, 5, 70, Appearance -> "Labeled"}
     ]

The vertices for that above can be generated by

    fbL = Flatten[Cases[Normal@FB, Line[x_] :> x, Infinity], 1]

Here is the original interface that I submitted, I will leave it in the question since there are several problems I had to solve in regard to manipulate and file saving that may be helpful to others.

    ClearAll;

    fx [i_, j_, d_, v_, h_, t_, w_, x0_, sf_] :=
          ( i (1 - d *Sin[v]) Cos[v] Sin[h] Sin[t] - 
          j (1 - d*Sin[v]) Cos[v] Cos[h] Cos[t] + 
          w*Log[2 - Sin[v]]*(-Sin[t] + Cos[t]))*sf + x0;


    fy [i_, j_, d_, v_, h_, t_, w_, y0_, sf_] :=
      (i (1 - d *Sin[v]) Cos[v] Sin[h] Cos[t] + 
          j (1 - d*Sin[v]) Cos[v] Cos[h] Sin[t] + 
          w*Log[2 - Sin[v]]*(Cos[t] + Sin[t]) )*sf + y0;


    fz[k_, v_, z0_, sf_] := 
      sf* k (1 - Sin[v])/2 + z0;


    fR[x_, y_, z_] := Sqrt[x^2 + y^2 + z^2];


    fL[x_, y_] := ArcTan[y/x];


    fB[x_, y_, z_] := ArcSin[z/Sqrt[x^2 + y^2 + z^2]];




    buildXYZCoords[np_, nq_, i_, j_, k_, d_, t_, w_, basex_] := 
     Module[{arr},
      pi = 3.14159265358979; 
      hstep = pi/(nq - 1); 
      vstep = 2*hstep;
      mb = fz[k, 0, 0, 1];
      sf = basex * Tan[d2r[mb]]/ mb;

      arr = {};
      For[v = 0, v <= 2 * pi, v = v + vstep,
       For[h =  -pi/2, h <= pi/2, h = h + hstep,
         x = fx[i, j, d, v, h, t, w, basex, sf];
         y = fy[i, j, d, v, h, t, w, 0, sf];
         z = fz[k, v, 0, sf];
         arr = Join[arr, {{x, y, z}}];
         ];
       ];
      (*  the last join must NOT have a semicolon, 
      a NULL return coccurs if it does, hence the line below *)
      arr = Join[arr, {}]
      ]

    buildCoords[np_, nq_, i_, j_, k_, d_, t_, w_, basex_] := Module[{arr},
      pi = N[Pi]; 
      hstep = pi/(nq - 1); 
      vstep = 2*hstep;
      mb = fz[k, 0, 0, 1];
      sf = basex * Tan[d2r[mb]]/ mb;

      arr = {{"id", "mid", "h", "v", "x", "y", "z", "R", "L", "B", "Ldeg",
          "Bdeg", "xm", "ym", "zm", "Rm", "Lm", "Bm", "Ldegm", "Bdegm"}};
      vid = 0;
      mid = 0;
      For[v = 0, v <= 2 * pi, v = v + vstep,
       For[h =  -pi/2, h <= pi/2, h = h + hstep,
         x = fx[i, j, d, v, h, t, w, basex, sf];
          y = fy[i, j, d, v, h, t, w, 0, sf];
          z = fz[k, v, 0, sf];
          R = fR[x, y, z];
          L = fL[x, y];
          B = fB[x, y, z];
          dL = r2d[L];
          dB = r2d[B];
          vid = vid + 1;

          arr = Join[arr, {{vid, mid, h, v, x, y, z, R, L, B, dL, dB}}];
         ];
       ];

      (* for export Semicolons as used above are required but the last \
    join must NOT have a semicolon i.e. 
      output not surpressed hence the line below *)
      arr = Join[arr, {}]
      ]


    r2d[r_] := N[r*180/Pi];

    d2r[d_] := N[d *Pi/180];



    m8 = Manipulate[
      FB = ParametricPlot3D[{
         fy [i, j, d, v, h, t, w, 0, 1],
         fx [i, j, d, v, h, t, w, k, 1],
         fz[k, v, 0, 1]},
        {v, 0, 2 π},
        {h, -π/2 , π/2},
        ViewPoint -> Dynamic@vp,
        Mesh -> Automatic,
        MeshStyle -> mcolor,
        PlotStyle -> scolor,
        AspectRatio -> 1,
        PlotRange -> {{-50, 50}, {0, 150}, {0, 100}},
        AxesLabel -> {"y", "x", "z"}
        ],


      "Parameters defining the surface shape",

      {{i, 35, "L ~y axis (i)"}, 5, 70, Appearance -> "Labeled"},
      {{j, 35, "To GC ~x axis (j)"}, 5, 70, Appearance -> "Labeled"},
      {{k, 50, "B ~z axis (k)" }, 5, 70, Appearance -> "Labeled"},
      {{d, 0.2, "Thinness (d)"}, 0, 1, Appearance -> "Labeled"},
      {{t, 0, "xy Tilt (t)"}, -Pi/2, Pi/2, Appearance -> "Labeled"},
      {{w, 0, "'Wind' factor (w)"}, -20, 20, Appearance -> "Labeled"},
      Delimiter,
      {{vp, {0, 8.5, 0}, "View Point"},
       {{8.5, 8.5, 8.5} -> "xyz", {0, 8.5, 0} -> "yz", {0, 0, 8.5} -> 
         "yx", {8.5, 0, 0} -> "xz"}, PopupMenu, ImageSize -> 140},
      Row[{Spacer[160], 
        Button["Reset", 
         FrontEndExecute[FrontEndToken["EvaluateNotebook"]], 
         ImageSize -> 140]}],
      {{mcolor, RGBColor[0.`, 1.`, 0.07`], "Mesh"}, ColorSlider},
      {{scolor, RGBColor[255, 72, 255], "Surface"}, ColorSlider},

      Delimiter,
      Row[
       {Spacer[100], Button["Export cartessian vertices", Export[
          SystemDialogInput["FileSave", 
           ToString[
            StringForm["FBdata_i=``_j=``_k=``_d=``_t=``_w=``.xlsx", i, j, 
             k, d, t, w]]], 
          TableForm[FB[[1, 1, 1]], 
           TableHeadings ->  {None, {"x", "y", "z"}}  ]]
         , Method -> "Queued", ImageSize -> 160]}],
      Row[{Spacer[100], Button["Export to STL", Export[
          SystemDialogInput["FileSave", 
           ToString[
            StringForm["FBdata_i=``_j=``_k=``_d=``_t=``_w=``.stl", i, j, 
             k, d, t, w]]], FB]
         , Method -> "Queued", ImageSize -> 160]}],

      Delimiter,
      "Export coordinates x,y,z,L,B,R with a specified density of data \
    points",
      "Distance of bubble base from observer",
      {{basex, 8.5, "Observer to Base: basex (Kpc)"}, 4, 10, 
       Appearance -> "Labeled"},
      {{ed, 60, "Export density"}, 10, 1000, Appearance -> "Labeled"},
      Row[{Spacer[100], 
        Button["Print Parameters", 
         Print["buildCoords[" <> ToString[ed] <> ","  <> ToString[ed] <> 
           "," <> ToString[i] <> "," <> ToString[j] <> "," <> 
           ToString[k] <> "," <> ToString[d] <> "," <> ToString[t] <> 
           "," <> ToString[w] <> "," <> ToString[basex] <> 
           "] // TableForm"],  Method -> "Queued", ImageSize -> 160]}],
      Row[{Spacer[100], 
        Button["Export Data", 
         Export[SystemDialogInput["FileSave", 
           ToString[
            StringForm[
             "FBdata_np=``_i=``_j=``_k=``_d=``_t=``_w=``_gcd=_``_full_\
    with_L_B.xlsx", ed, i, j, k, d, t, w, basex]]], 
          Evaluate@buildCoords[ed, ed, i, j, k, d, t, w, basex]], 
         Method -> "Queued", ImageSize -> 160]}],
      Row[{Spacer[100], 
        Button["Print Data", 
         Print[TableForm[
           Evaluate@buildCoords[ed, ed, i, j, k, d, t, w, basex]]], 
         Method -> "Queued", ImageSize -> 160]}],
      Row[{Spacer[100], 
        Button["Print Calculated XYZ", 
         Print[TableForm[
           Evaluate@buildXYZCoords[ed, ed, i, j, k, d, t, w, basex]]], 
         Method -> "Queued", ImageSize -> 160]}],
      Row[{Spacer[100], 
        Button["Plot Calculated Coords", 
         Print[ListSurfacePlot3D[
           buildXYZCoords[ed, ed, i, j, k, d, t, w, basex], 
           PlotRange -> {{0, 20}, {-10, 10}, {0, 20}}, 
           AxesLabel -> {"x", "y", "z"}, 
           PlotLabel -> {"Earth at (0,0,0)"}]], Method -> "Queued", 
         ImageSize -> 160]}]
      ]
$\endgroup$
  • 1
    $\begingroup$ You put really a lot of effort into the interface and I appreciate it. Alas, the actual problem could be presented with significantly less code. A true minimal example would make it much easier for others to isolate the problem and to come up with a tractable solution. $\endgroup$ – Henrik Schumacher Aug 5 at 13:20
  • $\begingroup$ I have simplified the question, I think that since Mathematica can display the surface it should not be difficult to work out intercepts, but no version of Nsolve of Findroot seems to be able to cope with the parametric equations $\endgroup$ – William Taylor Aug 9 at 9:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.