Return to mt_graph.rr CVS log

Up to [local] / OpenXM / src / asir-contrib / packages / src

Check the version to call pari(allocmem,10^8)

/* Two modules are defined:  mtg for plot3d and mtp for parametric_plot3d */ 
import("glib3.rr")$
if (version() < 20150801) {
  pari(allocatemem,10^8)$
}else{}$
/* pari(allocatemem,10^8)$
 BUG: 2015.08.07 it causes the cash of ox_pari server as a grandchild 
  when it is called in asirrc in ox_asir.
*/

module mtg;
static Tb$
static NN$
static Fmax$
static Xmin$
static Ymin$
static Xmax$
static Ymax$
static Trans$
static LevelMax$
static Slope$

static MZvalues$  /* value table of functions */
static Msize$     /* size of Mvmat */
static MXtoI$     /* (X-Xmin)*MXtoI --> index I */
static MYtoJ$     /* (Y-Ymin)*MYtoJ --> index J */
static MXfromI$
static MYfromJ$

static Gid$

Gid=0$

LevelMax = -1$  /* -1, no recursion */ 
Slope = 0.1$

#define MTSIZE 10

localf  plot3dInitTable;
localf  plot3d_nofit;
localf  plot3d;
localf  adaptive_outxyz;
localf  ueval;
localf  outxyz;
localf  outtexcoords;
localf  debug_mesh_init;
localf  debug_mesh;
localf  test1;
localf  test2;
localf  test3;

def plot3dInitTable() {
/*
  extern MZvalues$
  extern Msize$   
  extern MXtoI$   
  extern MYtoJ$   
  extern MXfromI$
  extern MYfromJ$

  extern LevelMax$
  extern NN$
  extern Xmax$
  extern Xmin$
  extern Ymax$
  extern Ymin$
*/

  if (LevelMax >= 0) {
    Msize = (NN+2)*(2^(LevelMax+1));
    MM = 2^(LevelMax+1);
  }
  else {
    Msize = (NN+2);
    MM = 1;
  }
  MZvalues=newmat(Msize,Msize);
  MXtoI = eval(exp(0)*NN*MM/(Xmax-Xmin));
  MYtoJ = eval(exp(0)*NN*MM/(Ymax-Ymin));
  MXfromI = 1/MXtoI;
  MYfromJ = 1/MYtoJ;
}

def plot3d_nofit(F) {
/*
  extern Tb;
  extern NN;
  extern Fmax;
  extern Trans;
*/

  Tb=string_to_tb("")$
  D=getopt(domain);
  if (type(D) <0) {
    D=[[-MTSIZE,MTSIZE],[-MTSIZE,MTSIZE]];
  }
  Xmin=D[0][0]; Xmax=D[0][1];
  Ymin=D[1][0]; Ymax=D[1][1];
  printf("domain=%a: Xmin=%a,Xmax=%a, Ymin=%a, Ymax=%a\n",D,Xmin,Xmax,Ymin,Ymax);
  if ((Xmin >= Xmax) || (Ymin >=Ymax)) error("Invalid domain.");
  /* (MTSIZE*2/(Xmax-Xmin))*(x-Xmin)-MTSIZE =a*x-b */
  /* (MTSIZE*2/(Ymax-Ymin))*(y-Ymin)-MTSIZE =c*y-d */
  Trans=newmat(2,2);  /* [[a,b],[c,d]] */
  Trans[0][0] = eval(exp(0)*MTSIZE*2/(Xmax-Xmin));
  Trans[0][1] = eval(-Xmin*exp(0)*MTSIZE*2/(Xmax-Xmin)-MTSIZE);
  Trans[1][0] = eval(exp(0)*MTSIZE*2/(Ymax-Ymin));
  Trans[1][1] = eval(-Ymin*exp(0)*MTSIZE*2/(Ymax-Ymin)-MTSIZE);
  /* print(Trans); */

  NN=getopt(mesh);
  if (type(NN) < 0) {
    NN=20;
  }
  printf("Mesh size is %a x %a. mesh=%a\n",NN,NN,NN);
  
  Xstep = eval(exp(0)*(Xmax-Xmin)/NN);
  Ystep = eval(exp(0)*(Ymax-Ymin)/NN);

  plot3dInitTable();


  /* Find the maximun of F */
  Fmax = ueval(F,Xmin,Ymin,-1,-1);  
  if (Fmax < 0) Fmax = -Fmax;
  Fmax_x = Xmin; Fmax_y = Ymin;
  for (X=Xmin; X<Xmax; X += Xstep) {
    for (Y=Ymin; Y<Ymax; Y += Ystep) {
      Z = ueval(F,X,Y,-1,-1);
      if (Z < 0)
        Z = -Z;
      if (Fmax < Z) {
        Fmax = Z;
        Fmax_x = X;
        Fmax_y = Y;
      }
    }
  }
  if (Fmax == 0.0) Fmax = 0.1;
  printf("Fmax is %a\n",Fmax);

  Vdollar = eval(exp(0)*(5-Trans[0][1])/Trans[0][0]); /* X coordinate */ 
  Vat = 5*Fmax/(MTSIZE*0.8);  /* Z coordinate.  outxyz() */ 
  printf("$=%a, @=%a (valid only for non-fit mode)\n",Vdollar,Vat);

  for (X=Xmin; X<Xmax; X += Xstep) {
    for (Y=Ymin; Y<Ymax; Y += Ystep) {
      adaptive_outxyz(F,X,Y,Xstep,Ystep,0,0);

    }
  }
  S=tb_to_string(Tb);

  Fname0 = "/tmp/oxm-"+getenv("USER")+"-";
  Fname = Fname0+rtostr(get_pid()+Gid)+".txt";
  Gid++;
  shell("rm -f "+Fname);
  output(Fname)$
  printf("%a",S)$
  output()$
  shell(getenv("OpenXM_HOME")+"/bin/oxmgraph "+Fname+" &");
  return [F,[["xmin",Xmin],["xmax",Xmax],["ymin",Ymin],["ymax",Ymax],
             ["mesh",NN],["fmax",Fmax]]];
}

def plot3d(F) {
/*
  extern Tb;
  extern NN;
  extern Fmax;
  extern Trans;
*/

  Fit=getopt(fit);
  if (type(Fit)!=1) Fit=0;
  if (Fit==0) return plot3d_nofit(F | option_list=getopt());

  /* Stupid of doing twice... should be rewritted. */
  D=getopt(domain);
  if (type(D) <0) {
    D=[[-MTSIZE,MTSIZE],[-MTSIZE,MTSIZE]];
  }
  Xmin=D[0][0]; Xmax=D[0][1];
  Ymin=D[1][0]; Ymax=D[1][1];
  if ((Xmin >= Xmax) || (Ymin >=Ymax)) error("Invalid domain.");

  NN=getopt(mesh);
  if (type(NN) < 0) {
    NN=20;
  }
  
  Xstep = eval(exp(0)*(Xmax-Xmin)/NN);
  Ystep = eval(exp(0)*(Ymax-Ymin)/NN);
  plot3dInitTable();

  /* Find the maximun of F */
  Fmax = ueval(F,Xmin,Ymin,-1,-1);  
  Vmax=Fmax; Vmin=Fmax;
  Fmax_x = Xmin; Fmax_y = Ymin; 
  Fmin_x = Xmin; Fmin_y = Ymin; 
  for (X=Xmin; X<Xmax; X += Xstep) {
    for (Y=Ymin; Y<Ymax; Y += Ystep) {
      Z = ueval(F,X,Y,-1,-1);
      if (Z > Vmax) { Vmax=Z; Fmax_x = X; Fmax_y = Y; }
      if (Z < Vmin) { Vmin=Z; Fmin_x = X; Fmin_y = Y; }
    }
  }
  Offset = -(Vmax+Vmin)/2;
  printf("Offset=%a\n",Offset);
  printf("Vmin=%a, [x,y]=%a\n",Vmin,[Fmin_x,Fmin_y]);
  printf("Vmax=%a, [x,y]=%a\n",Vmax,[Fmax_x,Fmax_y]);
  plot3d_nofit(F+Offset | option_list=getopt());
}

def adaptive_outxyz(F,X,Y,Xstep,Ystep,Xlevel,Ylevel) {
/*
  extern LevelMax;
  extern Slope;
  extern Fmax;
*/
  Xslope = (ueval(F,X,Y,-1,-1)-ueval(F,X+Xstep,Y,-1,-1))/Fmax; 
  if (Xslope < 0) Xslope = -Xslope;
  Yslope = (ueval(F,X,Y,-1,-1)-ueval(F,X,Y+Ystep,-1,-1))/Fmax; 
  if (Yslope < 0) Yslope = -Yslope;

  if ((Xlevel > LevelMax) || (Ylevel > LevelMax) || 
      ((Xslope < Slope) && (Yslope < Slope))) {
   Ystart=Y;
   for (I = 0; I < 2^Xlevel ; I++, X += Xstep) {
    for (J = 0, Y=Ystart; J < 2^Ylevel ; J++, Y += Ystep) {
      /* debug_mesh(X,Y,Xstep,Ystep); */
      write_to_tb("t ",Tb);      
      outxyz(F,X,Y); outxyz(F,X+Xstep,Y); outxyz(F,X,Y+Ystep);
      write_to_tb(" \n",Tb);
      write_to_tb("x ",Tb);
      outtexcoords(X,Y); outtexcoords(X+Xstep,Y); outtexcoords(X,Y+Ystep);
      write_to_tb(" \n",Tb);
      write_to_tb("t ",Tb);      
      outxyz(F,X,Y+Ystep); outxyz(F,X+Xstep,Y); outxyz(F,X+Xstep,Y+Ystep);
      write_to_tb(" \n",Tb);
      write_to_tb("x ",Tb);
      outtexcoords(X,Y+Ystep); outtexcoords(X+Xstep,Y); outtexcoords(X+Xstep,Y+Ystep);
      write_to_tb(" \n",Tb);
/*      Y += Ystep; */
    }

   }
  }else if ((Xslope >= Slope) && (Yslope >= Slope)) {
    adaptive_outxyz(F,X,Y,Xstep/2,Ystep/2,Xlevel+1,Ylevel+1);
  }else if ((Xslope < Slope) && (Yslope >= Slope)) {
    adaptive_outxyz(F,X,Y,Xstep,Ystep/2,Xlevel,Ylevel+1);
  }else if ((Xslope >= Slope) && (Yslope < Slope)) {
    adaptive_outxyz(F,X,Y,Xstep/2,Ystep,Xlevel+1,Ylevel);
  }else error("This case must not reached.");
}

def ueval(F,X,Y,I,J) {
/*
  extern MZvalues$
  extern Msize$   
  extern MXtoI$   
  extern MYtoJ$   
  extern MXfromI$
  extern MYfromJ$

  extern Xmax$
  extern Xmin$
  extern Ymax$
  extern Ymin$
*/

  /* Table lookup */
  if ((I >= 0) && (type(MZvalues[I][J]) != 0)) return MZvalues[I][J];

  if (I < 0) {
    I = pari(round,(X-Xmin)*MXtoI);
    J = pari(round,(Y-Ymin)*MYtoJ);
  } 
  if (type(MZvalues[I][J]) != 0) return MZvalues[I][J];

  if (type(F) == 17) { /* quote */
     Z = base_replace(F,[[x,X],[y,Y]]);
     Z = eval_quote(Z);
  }else{
     Z = subst(F,x,X,y,Y);
  }
  Z = eval(exp(0)*Z);
  MZvalues[I][J] = Z;
  return Z;
}
def outxyz(F,X,Y) {
/*
  extern Tb;
  extern Trans;
  extern Xmin;
  extern Ymin;
  extern MXtoI;
  extern MYtoJ;
  extern MXfromI;
  extern MYfromJ;
*/
  Z = ueval(F,X,Y,-1,-1);

  /* scaling */
  Z = Z*MTSIZE*0.8/Fmax;
  X = Trans[0][0]*X+Trans[0][1];
  Y = Trans[1][0]*Y+Trans[1][1];

  /* Truncate Z */
  if (Z > MTSIZE*2) Z = MTSIZE*2;
  if (Z < -MTSIZE*2) Z=-MTSIZE*2;

  /* Adjust numerical errors of X,Y */
  X = Xmin+pari(round,(X-Xmin)*MXtoI)*MXfromI;
  Y = Ymin+pari(round,(Y-Ymin)*MYtoJ)*MYfromJ;

  /* deval is used to avoid the expression like -8.1234 E-19 by Pari
  */
  S=sprintf(" %a %a %a ",deval(X),deval(Y),deval(Z));
  write_to_tb(S,Tb);
}

def outtexcoords(X,Y) {
/*
  extern Tb$
  extern NN$
  extern Trans$
  extern Xmin$
  extern Ymin$
  extern Xmax$
  extern Ymax$ 
*/
  N = NN/(Xmax-Xmin)*2;
  X = Trans[0][0]*X+Trans[0][1];
  Y = Trans[1][0]*Y+Trans[1][1];

  S = sprintf(" %a %a ",deval((X-Xmin)/(Xmax-Xmin)/N), deval((Y-Ymin)/(Ymax-Ymin)/N));
  write_to_tb(S,Tb);
}

def debug_mesh_init(X,Y) {
   glib_window(-X,-Y,X,Y);
   glib_clear();
}
/* Sample input. 
   test2();  
   LevelMax = 0;
   adaptive_outxyz(x,0.0, 0.0, 1.0, 1.0, 1,1);
*/
def debug_mesh(X,Y,Xstep,Ystep) {
   glib_line(X,Y,X+Xstep,Y);
   glib_line(X+Xstep,Y,X,Y+Ystep);
   glib_line(X,Y+Ystep,X,Y);

   glib_line(X,Y+Ystep,X+Xstep,Y);
   glib_line(X+Xstep,Y,X+Xstep,Y+Ystep);
   glib_line(X+Xstep,Y+Ystep,X,Y+Ystep);
   glib_flush();
} 

def test1() {
  plot3d((x^2+y^2)/30);
}

def test2() {
  /* debug_mesh_init(3,3); */
  plot3d(cos((x^2+y^2)^(1/2))+cos(3*(x^2+y^2)^(1/2)) | domain=[[-3,3],[-3,3]], mesh=30);
}

def test3() {
  plot3d(quote(sin(x)));
}
endmodule;


module mtp;
static Tb$
static NN$     /* mesh */
static Smin$
static Tmin$
static Smax$
static Tmax$

static Sxvalues$  /* value table of functions */
static Syvalues$  
static Szvalues$  
static Sevaluated$

static Pxmin$
static Pxmax$
static Pymin$
static Pymax$
static Pzmin$
static Pzmax$

#define MTSIZE 10    /* size of graph.c  [-MSIZE,MSIZE] */

localf parametric_plot3d_init;
localf parametric_plot3d;
localf p_outxyz;
localf myeval;
localf u3eval;
localf outxyz3;
localf test1;
localf test2;
localf test3;
localf test4;
localf test5;
localf test6;
localf test7;

def parametric_plot3d_init() {
/*
  extern NN$     

  extern Sxvalues$ 
  extern Syvalues$  
  extern Szvalues$  
  extern Sevaluated$
*/
  Sxvalues=newmat(NN+2,NN+2);
  Syvalues=newmat(NN+2,NN+2);
  Szvalues=newmat(NN+2,NN+2);
  Sevaluated=newmat(NN+2,NN+2);
}

def parametric_plot3d(F) {
/*
  extern Tb;
  extern NN;

  extern Smin$
  extern Tmin$
  extern Smax$
  extern Tmax$

  extern Pxmin$
  extern Pxmax$
  extern Pymin$
  extern Pymax$
  extern Pzmin$
  extern Pzmax$
*/

  Tb=string_to_tb("")$
  D=getopt(domain);
  if (type(D) <0) {
    D=[[-MTSIZE,MTSIZE],[-MTSIZE,MTSIZE]];
  }
  Smin=D[0][0]; Smax=D[0][1];
  Tmin=D[1][0]; Tmax=D[1][1];
  printf("domain=%a: Smin=%a,Smax=%a, Tmin=%a, Tmax=%a\n",D,Smin,Smax,Tmin,Tmax);
  if ((Smin >= Smax) || (Tmin >=Tmax)) error("Invalid domain.");

  NN=getopt(mesh);
  if (type(NN) < 0) {
    NN=20;
  }
  printf("Mesh size is %a x %a. mesh=%a\n",NN,NN,NN);
  
  Fitting=getopt(fitting);
  if (type(Fitting) < 0) {
    Fitting=1;
  }

  Sstep = eval(exp(0)*(Smax-Smin)/NN);
  Tstep = eval(exp(0)*(Tmax-Tmin)/NN);

  parametric_plot3d_init();

  /* Finding the bounding box */
  V = u3eval(F,Smin,Tmin);
  Pxmin = V[0]; Pxmax = V[0];
  Pymin = V[1]; Pymax = V[1];
  Pzmin = V[2]; Pzmax = V[2];
  for (S=Smin; S<Smax; S += Sstep) {
    for (T=Tmin; T<Tmax; T += Tstep) {
      V = u3eval(F,S,T);
      if (V[0] < Pxmin) Pxmin=V[0];
      if (V[1] < Pymin) Pymin=V[1];
      if (V[2] < Pzmin) Pzmin=V[2];
      if (V[0] > Pxmax) Pxmax=V[0];
      if (V[1] > Pymax) Pymax=V[1];
      if (V[2] > Pzmax) Pzmax=V[2];
    }
  }

  if (type(Fitting) == 0) {
    if (Pymin < Pxmin) Pxmin = Pymin;
    if (Pzmin < Pxmin) Pxmin = Pzmin; 
    if (Pymax > Pxmax) Pxmax = Pymax; 
    if (Pzmax > Pxmax) Pxmax = Pzmax;
    if (Pxmin < 0) Pxmin = - Pxmin;
    if (Pxmax < 0) Pxmax = - Pxmax;
    if (Pxmax < Pxmin) Pxmax = Pxmin;
    Pxmin = Pymin = Pzmin =  - Pxmax;
    Pxmax = Pymax = Pzmax =  Pxmax;
  }

  printf("xmin=%a, xmax=%a\n",Pxmin,Pxmax);
  printf("ymin=%a, ymax=%a\n",Pymin,Pymax);
  printf("zmin=%a, zmax=%a\n",Pzmin,Pzmax);

  /* cf. outxyz3() */
  Vdollar = (5+MTSIZE)*(Pxmax-Pymin)/(2*MTSIZE) + Pxmin;
  Vat =     (5+MTSIZE)*(Pzmax-Pzmin)/(2*MTSIZE) + Pzmin;
  printf("$=%a, @=%a\n",Vdollar,Vat);

  for (S=Smin; S<Smax; S += Sstep) {
    for (T=Tmin; T<Tmax; T += Tstep) {
      p_outxyz(F,S,T,Sstep,Tstep);
    }
  }
  SS=tb_to_string(Tb);
  Fname0 = "/tmp/oxm-"+getenv("USER")+"-";
  Fname = Fname0+rtostr(get_pid())+".txt";
  shell("rm -f "+Fname);
  output(Fname)$
  printf("%a",SS)$
  output()$
  shell(getenv("OpenXM_HOME")+"/bin/oxmgraph "+Fname+" &");
  return [F,[["xmin",Pxmin],["xmax",Pxmax],["ymin",Pymin],["ymax",Pymax],
             ["zmin",Pzmin],["zmax",Pzmax],
             ["mesh",NN],["dollar",Vdollar],["at",Vat]]];
}

def p_outxyz(F,S,T,Sstep,Tstep) { 
/*
  extern Tb$
  extern Smin$
  extern Tmin$
  extern Smax$
  extern Tmax$
*/
  write_to_tb("t ",Tb);      
  outxyz3(F,S,T); outxyz3(F,S+Sstep,T); outxyz3(F,S,T+Tstep);
  write_to_tb(" \n",Tb);

  write_to_tb("t ",Tb);      
  outxyz3(F,S,T+Tstep); outxyz3(F,S+Sstep,T); outxyz3(F,S+Sstep,T+Tstep);
  write_to_tb(" \n",Tb);
}

def myeval(Z) {
  return( eval(exp(0)*Z) ); 
}

def u3eval(F,S,T) {
/*
  extern NN$    
  extern Smin$
  extern Tmin$
  extern Smax$
  extern Tmax$

  extern Sxvalues$ 
  extern Syvalues$  
  extern Szvalues$  
  extern Sevaluated$
*/
  I = pari(round,(S-Smin)*NN/(Smax-Smin));
  J = pari(round,(T-Tmin)*NN/(Tmax-Tmin));

  if (type(Sevaluated[I][J]) != 0) 
   return [Sxvalues[I][J],Syvalues[I][J],Szvalues[I][J]];

  if (type(F[0]) == 17) { /* quote */
     Z = map(base_replace,F,[[s,S],[t,T]]);
     Z = map(eval_quote, Z);
  }else{
     Z = map(subst,F,s,S,t,T);
  }
  Z = map(myeval,Z);
  Sevaluated[I][J] = 1;
  Sxvalues[I][J] = Z[0]; 
  Syvalues[I][J] = Z[1]; 
  Szvalues[I][J] = Z[2]; 
  return Z;
}
def outxyz3(F,S,T) {
/*
  extern Tb;
  extern Smin;
  extern Tmin;
  extern NN;

  extern Pxmin$
  extern Pxmax$
  extern Pymin$
  extern Pymax$
  extern Pzmin$
  extern Pzmax$
*/
  V = u3eval(F,S,T);

  /* scaling */
  X = (V[0]-Pxmin)*MTSIZE*2/(Pxmax-Pxmin)-MTSIZE;
  Y = (V[1]-Pymin)*MTSIZE*2/(Pymax-Pymin)-MTSIZE;
  Z = (V[2]-Pzmin)*MTSIZE*2/(Pzmax-Pzmin)-MTSIZE;

  /* deval is used to avoid the expression like -8.1234 E-19 by Pari
  */
  S=sprintf(" %a %a %a ",deval(X),deval(Y),deval(Z));
  write_to_tb(S,Tb);
}


def test1() {
  parametric_plot3d([s,t,s^2+t^2]);
}

/* Sphere */
def test2() {
  parametric_plot3d([sin(s)*cos(t),sin(s)*sin(t),cos(s)] | domain=[[0,3.14],[0,3.14*2]]);
}

/* Mobius ring */
def test3() {
  parametric_plot3d([(2+t*cos(s/2))*cos(s),
                     (2+t*cos(s/2))*sin(s),
                     t*sin(s/2)] | domain=[[0,2*3.14],[-1,1]],
                                   fitting=0);
}

/* Torus */
def test4() {
  parametric_plot3d([(2+cos(t))*cos(s),
                     (2+cos(t))*sin(s),
                     sin(t)] | domain=[[0,2*3.14],[0,2*3.14]],
                               fitting=0);
}

/* Klein Bottle2,  Wikipedia, Figure 8 immersion*/
def test5() {
  L=    [5*cos(x)*(cos(x/2)*cos(y)+sin(x/2)*sin(2*y)+3.0)-10.0,
         -5*sin(x)*(cos(x/2)*cos(y)+sin(x/2)*sin(2*y)+3.0),
         5*(-sin(x/2)*cos(y)+cos(x/2)*sin(2*y))];
  L = map(subst,L,x,s,y,t);
  parametric_plot3d(L | domain=[[-3.14,3.14],[-3.14,3.14]], fitting=0);
}

/* Klein Bottle, k3dsurf 'Klein demo' */
def test6() {
  L=    [(3*(1+sin(y)) + 2*(1-cos(y)/2)*cos(x))*cos(y),
	/*cos(x)*(cos(x/2)*(2^(1/2)+cos(y))+sin(x/2)*sin(y)*cos(y)),*/
	(4+2*(1-cos(y)/2)*cos(x))*sin(y),
	/*sin(x)*(cos(x/2)*(2^(1/2)+cos(y))+sin(x/2)*sin(y)*cos(y)),*/
	-2*(1-cos(y)/2) * sin(x)
	/*(-1)*sin(x/2)*(2^(1/2)+cos(y))+cos(x/2)*sin(y)*cos(y)*/
	];
  L = map(subst,L,x,s,y,t);
  parametric_plot3d(L | domain=[[0.0,2*3.14],[0.0,2*3.14]]);
}

def test7() {
   parametric_plot3d([s,s^2+t/10, s^3+t/10]);
}

endmodule;

end$