5

I have a surface specified by samples of a function on a square grid (in particular, there is no formula for the function), and I want to draw a surface plot along with some contour lines.

As you can see in the image below, just calculating contour lines with the contour package doesn't work because the contour calculation doesn't account for the surface plot smoothing. (Actually, it doesn't seem to work even without smoothing, i.e., dropping Spline and using -- as the join operator for the contour lines).

import graph3;
import contour; 
import palette;
import x11colors; 
size3(200,IgnoreAspect);
settings.outformat="pdf";
settings.render=8;
currentprojection=perspective(20,8,2); 

// A function specified by its values on a 10 x 10 grid
real[][] z =
  {{0,0,0,0,0,0,0,0,0,0},
   {0,0.459792676157358,0.136148250882153,-0.163615908155485,-0.0259463996973455,-0.0840769425808095,0.62126438834251,0.518640203087028,1.03569619721556,0},
   {0,-0.604162941996047,-0.424941391270887,0.259335710132218,0.650164067086578,0.127505396459818,0.302247026415486,0.254979306532084,0.252221600145993,0},
   {0,0.0304849331645697,-0.524980730553991,0.854096681304495,0.480467910001354,0.258129875403682,0.692548824372049,-0.252410071666853,-0.167082932245632,0},
   {0,0.987911902753579,0.382866021557655,0.864790291964691,0.607114025159436,0.55793941674487,0.441528479985084,0.767320637646196,0.467364706957709,0},
   {0,0.115108721902316,0.516912837384094,1.01635622832271,0.660704188001264,0.382332982489028,-0.307456859893877,-0.382304124975139,1.06361557239264,0},
   {0,0.288185323850592,-0.117203025712584,0.803252469880009,0.714506160677705,0.151669905008431,-0.342341661772393,-0.0261977699813443,0.473208218671247,0},
   {0,0.726871540197522,0.0683933588314665,0.481762795614408,0.969192479689294,-0.40613837131098,0.142722121638684,0.635132993615454,0.237704962254851,0},
   {0,-0.769077284482047,-0.755640443336741,0.892723307000915,0.683470361768723,0.470367719595602,0.967570149005728,0.249827942333624,0.478612319838185,0},
   {0,0,0,0,0,0,0,0,0,0}};
// Dimensions of the array
int m = 10, n = 10; 

triple f(pair t) {
  int i=round(t.x);
  int j=round(t.y);
  return (i,j,z[i][j]);
}

// color and draw surface
surface s=surface(f,(0,0),(m-1,n-1),m-1,n-1,Spline);
s.colors(palette(s.map(zpart),Rainbow()));
draw(s); 

// draw a contour line with height 0.25 
real t = 0.25; 
guide[][] contours = contour(z,(0,0),(9,9),new real[]{t},join=operator ..);
guide3[][] liftedcontours = lift(new real (real x, real y) {return t;}, contours);
draw(liftedcontours[0][1],linewidth(2.0));

Surface plot with contour

6

z is only discrete and there is no reason that the (smooth) contours lines of z are the contours lines of the smooth surface. You have to construct the function corresponding to the surface, it is possible through s.point and then you can use contour. Please find a possible code

import graph3;
import contour; 
import palette;
import x11colors; 
size3(200,IgnoreAspect);
//settings.outformat="pdf";
settings.render=8;
currentprojection=perspective(20,8,2); 

// A function specified by its values on a 10 x 10 grid
real[][] z =
  {{0,0,0,0,0,0,0,0,0,0},
   {0,0.459792676157358,0.136148250882153,-0.163615908155485,-0.0259463996973455,-0.0840769425808095,0.62126438834251,0.518640203087028,1.03569619721556,0},
   {0,-0.604162941996047,-0.424941391270887,0.259335710132218,0.650164067086578,0.127505396459818,0.302247026415486,0.254979306532084,0.252221600145993,0},
   {0,0.0304849331645697,-0.524980730553991,0.854096681304495,0.480467910001354,0.258129875403682,0.692548824372049,-0.252410071666853,-0.167082932245632,0},
   {0,0.987911902753579,0.382866021557655,0.864790291964691,0.607114025159436,0.55793941674487,0.441528479985084,0.767320637646196,0.467364706957709,0},
   {0,0.115108721902316,0.516912837384094,1.01635622832271,0.660704188001264,0.382332982489028,-0.307456859893877,-0.382304124975139,1.06361557239264,0},
   {0,0.288185323850592,-0.117203025712584,0.803252469880009,0.714506160677705,0.151669905008431,-0.342341661772393,-0.0261977699813443,0.473208218671247,0},
   {0,0.726871540197522,0.0683933588314665,0.481762795614408,0.969192479689294,-0.40613837131098,0.142722121638684,0.635132993615454,0.237704962254851,0},
   {0,-0.769077284482047,-0.755640443336741,0.892723307000915,0.683470361768723,0.470367719595602,0.967570149005728,0.249827942333624,0.478612319838185,0},
   {0,0,0,0,0,0,0,0,0,0}};
// Dimensions of the array
int m = 10, n = 10; 

triple f(pair t) {
  int i=round(t.x);
  int j=round(t.y);
  return (i,j,z[i][j]);
}

// color and draw surface
surface s=surface(f,(0,0),(m-1,n-1),m-1,n-1,Spline);
s.colors(palette(s.map(zpart),Rainbow()));
draw(s);
real g(real a, real b)
{
return s.point(a,b).z;
}
// draw a contour line with height 0.25 
real t = 0.25; 
guide[][] contours = contour(g,(0,0),(8.9,8.9),new real[]{t},50, join=operator ..);
//guide[][] contours = contour(z,(0,0),(9,9),new real[]{t},join=operator ..);
guide3[][] liftedcontours = lift(new real (real x, real y) {return t;}, contours);
draw(liftedcontours[0][1],linewidth(2.0));

And the pictureenter image description here

Notice that in contour(g,(0,0),(8.9,8.9),new real[]{t},50, join=operator ..); if 8.9 is replaced by 9 then an index error occurs. 50 (nx=50) is to reduce the accuracy, it is sufficient. A lower value produces also an index error...

O.G.

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.