Asymptote using graph.asy – fig0230

Category: Asymptote,Examples 2D,graph.asyPh. Ivaldi @ 0 h 10 min

Figure 0022
(Compiled with Asymptote version 2.14svn-r5318)
    
import graph;
unitsize(x=1cm,y=2cm);

struct rational
{
  int p;
  int q;
  real ep=1/10^5;
};

rational operator init() {return new rational;}

rational rational(real x, real ep=1/10^5)
{
  rational orat;
  int q=1;
  while (abs(round(q*x)-q*x)>ep)
    {
      ++q;
    }
  orat.p=round(q*x);
  orat.q=q;
  orat.ep=ep;
  return orat;
}

int pgcd(int a, int b)
{
  int a_=abs(a), b_=abs(b), r=a_;
  if (b_>a_) {a_=b_; b_=r; r=a_;}
  while (r>0)
    {
      r=a_%b_;
      a_=b_;
      b_=r;
    }
  return a_;
}

string texfrac(int p, int q,
               string factor="",
               bool signin=false, bool factorin=true,
               bool displaystyle=false,
               bool zero=true)
{
  if (p==0) return (zero ? "$0$" : "");
  string disp= displaystyle ? "$\displaystyle " : "$";
  int pgcd=pgcd(p,q);
  int num= round(p/pgcd), den= round(q/pgcd);
  string nums;
  if (num==1)
    if (factor=="" || (!factorin && (den !=1))) nums="1"; else nums="";
  else
    if (num==-1)
      if (factor=="" || (!factorin && (den !=1))) nums="-1"; else nums="-";
    else nums= (string) num;
  if (den==1) return "$" + nums + factor + "$";
  else
    {
      string dens= (den==1) ? "" : (string) den;
      if (signin || num>0)
        if (factorin)
          return disp + "\frac{" + nums + factor + "}{" + (string) dens + "}$";
        else
          return disp + "\frac{" + nums + "}{" + (string) dens + "}"+ factor + "$";
      else
        {
          if (num==-1)
            if (factor=="" || !factorin) nums="1"; else nums="";
          else nums=(string)(abs(num));
        if (factorin)
          return disp + "-\frac{" + nums + factor + "}{" + (string) dens + "}$";
        else
          return disp + "-\frac{" + nums + "}{" + (string) dens + "}"+ factor + "$";
        }
    }
}

string texfrac(rational x,
               string factor="",
               bool signin=false, bool factorin=true,
               bool displaystyle=false,
               bool zero=true)
{
  return texfrac(x.p, x.q, factor, signin, factorin, displaystyle, zero);
}

ticklabel labelfrac(real ep=1/10^5, real factor=1.0,
                    string symbol="",
                    bool signin=false, bool symbolin=true,
                    bool displaystyle=false,
                    bool zero=true)
{
  return new string(real x)
    {
      return texfrac(rational(x/factor), symbol, signin, symbolin, displaystyle, zero);
    };
}

ticklabel labelfrac=labelfrac();

xlimits( -2pi, 2pi);
ylimits( -1, 1);

yaxis("y",LeftRight , Ticks(labelfrac,Step=.5,step=.25, ptick=grey, extend=true));

xaxis("$\theta$",BottomTop, Ticks(labelfrac(factor=pi,symbol="\pi",symbolin=false),
                           Step=pi/2, step=pi/4, ptick=grey, extend=true));

draw(graph(new real(real x){return sin(x);},-2pi,2pi));
draw(graph(new real(real x){return cos(x);},-2pi,2pi), .8red);

Étiquettes : , , , , ,


Asymptote using graph.asy – fig0240

Category: Asymptote,Examples 2D,graph.asyPh. Ivaldi @ 1 h 10 min

Figure 0023
(Compiled with Asymptote version 1.84svn-r4619)
    
import graph;

// public real xunit=1cm,yunit=1cm;

void graphicrules(picture pic=currentpicture, string prefix=defaultfilename, real unit=1cm,
                  real xunit=unit != 0 ? unit : 0,
                  real yunit=unit != 0 ? unit : 0,
                  real xmin, real xmax, real ymin, real ymax)
{
  xlimits(xmin, xmax);
  ylimits(ymin, ymax);
  unitsize(x=xunit, y=yunit);
}

struct rational
{
  int p;
  int q;
  real ep=1/10^5;
};

rational operator init() {return new rational;}

rational rational(real x, real ep=1/10^5)
{
  rational orat;
  int q=1;
  while (abs(round(q*x)-q*x)>ep)
    {
      ++q;
    }
  orat.p=round(q*x);
  orat.q=q;
  orat.ep=ep;
  return orat;
}

int pgcd(int a, int b)
{
  int a_=abs(a), b_=abs(b), r=a_;
  if (b_>a_) {a_=b_; b_=r; r=a_;}
  while (r>0)
    {
      r=a_%b_;
      a_=b_;
      b_=r;
    }
  return a_;
}

string texfrac(int p, int q,
               string factor="",
               bool signin=false, bool factorin=true,
               bool displaystyle=false,
               bool zero=true)
{
  if (p==0) return (zero ? "$0$" : "");
  string disp= displaystyle ? "$\displaystyle " : "$";
  int pgcd=pgcd(p,q);
  int num= round(p/pgcd), den= round(q/pgcd);
  string nums;
  if (num==1)
    if (factor=="" || (!factorin && (den !=1))) nums="1"; else nums="";
  else
    if (num==-1)
      if (factor=="" || (!factorin && (den !=1))) nums="-1"; else nums="-";
    else nums= (string) num;
  if (den==1) return "$" + nums + factor + "$";
  else
    {
      string dens= (den==1) ? "" : (string) den;
      if (signin || num>0)
        if (factorin)
          return disp + "\frac{" + nums + factor + "}{" + (string) dens + "}$";
        else
          return disp + "\frac{" + nums + "}{" + (string) dens + "}"+ factor + "$";
      else
        {
          if (num==-1)
            if (factor=="" || !factorin) nums="1"; else nums="";
          else nums=(string)(abs(num));
          if (factorin)
            return disp + "-\frac{" + nums + factor + "}{" + (string) dens + "}$";
          else
            return disp + "-\frac{" + nums + "}{" + (string) dens + "}"+ factor + "$";
        }
    }
}

string texfrac(rational x,
               string factor="",
               bool signin=false, bool factorin=true,
               bool displaystyle=false,
               bool zero=true)
{
  return texfrac(x.p, x.q, factor, signin, factorin, displaystyle, zero);
}

ticklabel labelfrac(real ep=1/10^5, real factor=1.0,
                    string symbol="",
                    bool signin=false, bool symbolin=true,
                    bool displaystyle=false,
                    bool zero=true)
{
  return new string(real x)
    {
      return texfrac(rational(x/factor), symbol, signin, symbolin, displaystyle, zero);
    };
}

ticklabel labelfrac=labelfrac();

void grid(picture pic=currentpicture,
          real xmin=pic.userMin().x, real xmax=pic.userMax().x,
          real ymin=pic.userMin().y, real ymax=pic.userMax().y,
          real xStep=1, real xstep=.5,
          real yStep=1, real ystep=.5,
          pen pTick=nullpen, pen ptick=grey, bool above=true)
{
  xaxis(pic, BottomTop, xmin, xmax, Ticks("%",extend=true,Step=xStep,step=xstep,pTick=pTick,ptick=ptick), above=above);
  yaxis(pic, LeftRight, ymin, ymax, Ticks("%",extend=true,Step=yStep,step=ystep,pTick=pTick,ptick=ptick), above=above);
}

void cartesianaxis(picture pic=currentpicture,
                   Label Lx=Label("$x$",align=S),
                   Label Ly=Label("$y$",align=W),
                   real xmin=pic.userMin().x, real xmax=pic.userMax().x,
                   real ymin=pic.userMin().y, real ymax=pic.userMax().y,
                   real extrawidth=1, real extraheight=extrawidth,
                   pen p=currentpen,
                   ticks xticks=Ticks("%",pTick=nullpen, ptick=grey),
                   ticks yticks=Ticks("%",pTick=nullpen, ptick=grey),
                   bool above=true,
                   arrowbar arrow=Arrow)
{
  extraheight= cm*extraheight/(2*pic.yunitsize);
  extrawidth = cm*extrawidth/(2*pic.xunitsize);
  yequals(pic, Lx, 0, xmin-extrawidth, xmax+extrawidth, p, above, arrow=arrow);
  yequals(pic, 0, xmin, xmax, p, xticks, above);
  xequals(pic, Ly, 0, ymin-extraheight, ymax+extraheight, p, above, arrow=arrow);
  xequals(pic, 0, ymin, ymax, p, yticks, above);
}

void labeloij(picture pic=currentpicture,
              Label Lo=Label("$O$",NoFill),
              Label Li=Label("$\vec{\imath}$",NoFill),
              Label Lj=Label("$\vec{\jmath}$",NoFill),
              pair diro=SW, pair diri=S, pair dirj=W,
              pen p=scale(2)*currentpen,
              filltype filltype=NoFill, arrowbar arrow=Arrow(2mm))
{
  if (Lo.filltype==NoFill) Lo.filltype=filltype;
  if (Li.filltype==NoFill) Li.filltype=filltype;
  if (Lj.filltype==NoFill) Lj.filltype=filltype;
  labelx(pic, Lo, 0, diro, p);
  draw(pic, Li, (0,0)--(1,0), diri, p, arrow);
  draw(pic, Lj, (0,0)--(0,1), dirj, p, arrow);
  dot(pic, (0,0), dotsize(p)+p);
}

void labeloIJ(picture pic=currentpicture,
              Label Lo=Label("$O$",NoFill),
              Label LI=Label("$I$",NoFill),
              Label LJ=Label("$J$",NoFill),
              pair diro=SW, pair dirI=S, pair dirJ=W,
              pen p=currentpen,
              filltype filltype=NoFill, arrowbar arrow=Arrow)
{
  if (Lo.filltype==NoFill) Lo.filltype=filltype;
  if (LI.filltype==NoFill) LI.filltype=filltype;
  if (LJ.filltype==NoFill) LJ.filltype=filltype;
  labelx(pic, LI, 1, dirI, p);
  labely(pic, LJ, 1, dirJ, p);
  labelx(pic, Lo, 0, diro, p);
  dot(pic, (0,0), dotsize(p)+p);
}

graphicrules(xunit=1cm, yunit=3cm,
             xmin=-2pi, xmax=2pi, ymin=-1, ymax=1);
grid(xStep=pi/2, xstep=pi/4, yStep=.5, ystep=.25);
cartesianaxis(xticks=Ticks(Label(UnFill),labelfrac(factor=pi,symbol="\pi",symbolin=true, zero=false),Step=pi/2, step=pi/4, ptick=grey),
              yticks=Ticks(Label(UnFill),labelfrac(zero=false),Step=.5,step=.25, ptick=grey), arrow=None);
dot("$O$",(0,0),2SW);

Étiquettes : , , , , ,


Random walk in the space – fig0020

Category: Asymptote,Examples 3D,Random Walk (3D)Ph. Ivaldi @ 12 h 46 min

Figure 0002
(Compiled with Asymptote version 1.84svn-r4619)
    
import three;
settings.render=0;

// The available directions of steps
triple[] dirs={X,-X,Y,-Y,Z,-Z};
dirs.cyclic=true;

struct walk
{
  triple[] nodes;
  pen[] p;
}

// Comput the nodes of the path
walk randWalk(real Srnd(), int n, pen[] p={currentpen})
{
  p.cyclic=true;
  walk ow;
  triple pos=O, tpos;
  for (int i=0; i < n; ++i) {
    int R=round(Srnd());
    tpos=pos+dirs[R];
    ow.nodes.push(tpos);
    ow.p.push(p[R]);
    pos=tpos;
  }
  return ow;
}

walk randWalk(int Srnd(), int n, pen[] p={currentpen})
{
  real R(){ return Srnd();}
  return randWalk(R,n,p);
}

void drawWalk(walk walk)
{
  triple camera=currentprojection.camera;
  if(currentprojection.infinity)
    camera *= max(abs(minbound(walk.nodes)),abs(maxbound(walk.nodes)));
  real[][] depth;
  for(int i=0; i < walk.nodes.length-1; ++i) {
    real d=abs(camera-0.5*(walk.nodes[i]+walk.nodes[i+1]));
    depth.push(new real[] {d,i});
  }
  depth=sort(depth);
  triple M=walk.nodes[round(depth[0][1])];
  triple m=walk.nodes[round(depth[depth.length-1][1]+1)];
  // Draw from farthest to nearest
  while(depth.length > 0) {
    real[] a=depth.pop();
    int i=round(a[1]);
    // dot(walk.nodes[i],walk.p[i]);
    draw(walk.nodes[i]--walk.nodes[i+1],abs(walk.nodes[i]-m)/abs(M-m)*(walk.p[i]+walk.p[i+1]));
  }
}


size(18cm);
currentprojection=orthographic((0.5,0.5,1));

drawWalk(randWalk(rand,50000,new pen[]{red, blue, green, yellow, purple}));
shipout(bbox(3mm,Fill));

Étiquettes : ,


Official Asymptote example – elliptic

Category: Asymptote,Official Gallery One-PagerPh. Ivaldi @ 18 h 57 min

Figure 0054
(Compiled with Asymptote version 2.14svn-r5318)
/* This code comes from The Official Asymptote Gallery */
    
struct curve {
  real a=0;
  real b=8;

  real y2(real x) {
    return x^3+a*x+b;
  }

  real disc() {
    return -16*(4*a*a*a+27*b*b);
  }

  real lowx () {
    return sqrt(-a/3);
  }

  int comps() {
    if (a < 0) {
      real x=sqrt(-a/3);
      return y2(x) < 0 ? 2 : 1;
    }
    return 1;
  }

  void locus(picture pic=currentpicture, real m, real M, int n=100,
             pen p=currentpen) {
    path flip(path p, bool close) {
      path pp=reverse(yscale(-1)*p)..p;
      return close ? pp..cycle : pp;
    }
    path section(real m, real M, int n) {
      guide g;
      real width=(M-m)/n;
      for(int i=0; i <= n; ++i) {
        real x=m+width*i;
        real yy=y2(x);
        if (yy > 0)
          g=g..(x,sqrt(yy));
      }
      return g;
    }

    if (comps() == 1) {
      draw(pic,flip(section(m,M,n),false),p);
    }
    else {
      real x=lowx(); // The minimum on x^3+ax+b
      if (m < x)
        draw(pic,flip(section(m,min(x,M),n),true),p);
      if (x < M)
        draw(pic,flip(section(max(x,m),M,n),false),p);
    }
  }

  pair neg(pair P) {
    return finite(P.y) ? yscale(-1)*P : P;
  }

  pair add(pair P, pair Q) {
    if (P.x == Q.x && P.x != Q.x)
      return (0,infinity);
    else {
      real lambda=P == Q ? (3*P.x^2+a)/(2*P.y) : (Q.y-P.y)/(Q.x-P.x);
      real Rx=lambda^2-P.x-Q.x;
      return (Rx,(P.x-Rx)*lambda-P.y);
    }
  }
}

import graph;
import math;

size(0,200);

curve c; c.a=-1; c.b=4;

pair oncurve(real x) 
{
  return (x,sqrt(c.y2(x)));
}

picture output;

axes();
c.locus(-4,3,.3red+.7blue);

pair P=oncurve(-1),Q=oncurve(1.2);
pair PP=c.add(P,P),sum=c.add(P,Q);

save();

drawline(P,Q,dashed);
drawline(c.neg(sum),sum,dashed);
dot("$P$", P, NW);
dot("$Q$", Q, SSE);
dot(c.neg(sum));
dot("$P+Q$", sum, 2SW);

add(output,currentpicture.fit(),(-0.5cm,0),W);

restore();

save();

drawline(P,c.neg(PP),dashed);
drawline(c.neg(PP),PP,dashed);
dot("$P$", P, NW);
dot(c.neg(PP));
dot("$2P$", PP, SW);

add(output,currentpicture.fit(),(0.5cm,0),E);

shipout(output);
    
restore();

Étiquettes : , ,


Official Asymptote example – pdb

Category: Asymptote,Official Gallery One-PagerPh. Ivaldi @ 14 h 57 min

Figure 0152
(Compiled with Asymptote version 1.87svn-r4652)
/* This code comes from The Official Asymptote Gallery */
    
import three;
import cpkcolors;

// A sample Protein Data Bank file for this example is available from
// http://ndbserver.rutgers.edu/ftp/NDB/coordinates/na-biol/100d.pdb1

bool getviews=true;

currentlight=White;
//currentlight=nolight;

size(200);
currentprojection=perspective(30,30,15);

// Uncomment this line for more accurate (but slower) PDF rendering
//dotgranularity=0;

pen chainpen=green;
pen hetpen=purple;

string filename="100d.pdb1";
//string filename=getstring("filename");

string prefix=stripextension(filename);
file data=input(filename);

pen color(string e) 
{
  e=replace(e," ","");
  int n=length(e);
  if(n < 1) return currentpen;
  if(n > 1) e=substr(e,0,1)+downcase(substr(e,1,n-1));
  int index=find(Element == e);
  if(index < 0) return currentpen;
  return rgb(Hexcolor[index]);
}       

// ATOM
string[] name,altLoc,resName,chainID,iCode,element,charge;
int[] serial,resSeq;
real[][] occupancy,tempFactor;

bool newchain=true;

struct bond 
{
  int i,j;
  void operator init(int i, int j) {
    this.i=i;
    this.j=j;
  }
}

bond[] bonds;

struct atom 
{
  string name;
  triple v;
  void operator init(string name, triple v) {
    this.name=name;
    this.v=v;
  }
}

struct chain
{
  int[] serial;
  atom[] a;
}

int[] serials;
chain[] chains;
atom[] atoms;

while(true) {
  string line=data;
  if(eof(data)) break;
  string record=replace(substr(line,0,6)," ","");
  if(record == "TER") {newchain=true; continue;}
  bool ATOM=record == "ATOM";
  bool HETATOM=record == "HETATM";
  int serial;

  atom a;
  if(ATOM || HETATOM) {
    serial=(int) substr(line,6,5);
    a.name=substr(line,76,2);
    a.v=((real) substr(line,30,8),
         (real) substr(line,38,8),
         (real) substr(line,46,8));
  }
  if(ATOM) {
    if(newchain) {
      chains.push(new chain);
      newchain=false;
    }
    chain c=chains[chains.length-1];
    c.serial.push(serial);
    c.a.push(a);
    continue;
  }
  if(HETATOM) {
    serials.push(serial);
    atoms.push(a);
  }
  if(record == "CONECT") {
    int k=0;
    int i=(int) substr(line,6,5);
    while(true) {
      string s=replace(substr(line,11+k,5)," ","");
     if(s == "") break;
      k += 5;
      int j=(int) s;
      if(j <= i) continue;
      bonds.push(bond(i,j));
     }
  }
}

write("Number of atomic chains: ",chains.length);

int natoms;
for(chain c : chains) {
  for(int i=0; i < c.a.length-1; ++i)
    draw(c.a[i].v--c.a[i+1].v,chainpen,currentlight);
  for(atom a : c.a)
    dot(a.v,color(a.name),currentlight);
  natoms += c.a.length;
}

write("Number of chained atoms: ",natoms);
write("Number of hetero atoms: ",atoms.length);

for(atom h : atoms)
  dot(h.v,color(h.name),currentlight);

write("Number of hetero bonds: ",bonds.length);

for(bond b : bonds) {
  triple v(int i) {return atoms[find(serials == i)].v;}
  draw(v(b.i)--v(b.j),hetpen,currentlight);
}

string options;
string viewfilename=prefix+".views";

if(!error(input(viewfilename,check=false)))
  options="3Dviews="+viewfilename;

if(getviews) {
  picture pic;
  add(pic,embed("label",currentpicture,options=options),(0,0),N);
  label(pic,cameralink("label"),(0,0),S,fontsize(12pt));
  shipout(prefix,pic,options=options);
} else
  shipout(prefix,options=options);

Étiquettes : , , , , ,


Various code with Asymptote – fig0300

Category: Asymptote,MiscellaneousPh. Ivaldi @ 23 h 11 min

Figure 0003
(Compiled with Asymptote version 1.87svn-r4652)
    
// On définit ce qu'est l'objet 'rectangle'
// On peut comparer ça à la définition d'une voiture:
// elle a des propriètés: la couleur, le nombre de portes...
// elle a aussi des actions (ou méthodes en info): ouvrir la porte, tourner le volant, avancer...
struct rectangle
{
  pair A,B,C,D; // Ce sont des propriétés (au sens informatique): ses sommets.
  path path(){return A--B--C--D--cycle;} // Ceci est une action du rectangle (renvoyer le chemin du rectangle)
};

// Le constructeur de l'objet (une action qui renvoie l'objet).
// Chaque fois que l'on va définir une variable 'X' de type 'rectangle'
// c'est le retour de l'action rectangle.init qui est affecté à X.
rectangle operator init() {return new rectangle;}

// La fonction rectangle  qui prend en argument le coin inférieur gauche
// et le coin supérieur droit retourne un rectangle.
// On notera comment l'on récupére une propriété ou une action d'un
// objet grâce au '.' en écrivant objet.propriété ou
// objet.action(les_paramètres).
rectangle rectangle(pair gb, pair dh)
{
  rectangle retour;
  retour.A=(dh.x,gb.y);
  retour.B=dh;
  retour.C=(gb.x,dh.y);
  retour.D=gb;
  return retour;
}

// On peut alors définir tous les opérateurs que l'on veut sur les
// rectangles.
// par exemple la multiplication d'un réal par un rectangle (x*rectangle):
rectangle operator *(real x, rectangle rectangle)
{
  rectangle retour;
  retour.A=rectangle.A+x*(rectangle.B-rectangle.A);
  retour.B=rectangle.B+x*(rectangle.C-rectangle.B);
  retour.C=rectangle.C+x*(rectangle.D-rectangle.C);
  retour.D=rectangle.D+x*(rectangle.A-rectangle.D);
  return retour;
}

// Ici, cette opération est commutative (rectangle*x=x*rectangle):
rectangle operator *(rectangle rectangle,real x){return x*rectangle;}

// On dit ici comment "Asymptote" doit convertir, si besoin est, un type 'rectangle'
// en 'path' puis en 'path[]' (tableau de 'path').
path operator cast(rectangle rectangle){return rectangle.path();}
path[] operator cast(rectangle rectangle){return rectangle.path();}
//Par exemple, dans le code:
/*
  rectangle SQ=rectangle((0,0), (1,1));
  draw(SQ);
*/
// 'draw' attend un type 'path' donc 'SQ' sera automatiquement converti.


// Maintenant on peut travailler avec cette structure
size(10cm);
rectangle SQ=rectangle((0,0), (1,1));
pen [] col= new pen[]{gray,yellow};

int nb=10;
for (int i=0; i<nb; ++i) {
  filldraw(SQ,col[i%2]);
  SQ=.25*SQ;
 }

Étiquettes :