Asymptote using graph3.asy – fig0180

Category: Asymptote,Examples 3D,graph3.asyPh. Ivaldi @ 20 h 11 min

Figure 0018

The spherical harmonics latex2png equation are the angular portion of the solution to Laplace's equation in spherical coordinates where azimuthal symmetry is not present.

The spherical harmonics are defined by:

latex2png equation

where latex2png equation and latex2png equation is the Legendre polynomial.

Source

(Compiled with Asymptote version 2.14svn-r5318)
    
import palette;
import math;
import graph3;

typedef real fct(real);
typedef pair zfct2(real,real);
typedef real fct2(real,real);

real binomial(real n, real k)
{
  return gamma(n+1)/(gamma(n-k+1)*gamma(k+1));
}

real factorial(real n) {
  return gamma(n+1);
}

real[] pdiff(real[] p)
{ // p(x)=p[0]+p[1]x+...p[n]x^n
  // retourne la dérivée de p
  real[] dif;
  for (int i : p.keys) {
    if(i != 0) dif.push(i*p[i]);
  }
  return dif;
}

real[] pdiff(real[] p, int n)
{ // p(x)=p[0]+p[1]x+...p[n]x^n
  // dérivée n-ième de p
  real[] dif={0};
  if(n >= p.length) return dif;
  dif=p;
  for (int i=0; i < n; ++i)
    dif=pdiff(dif);
  return dif;
}

fct operator *(real y, fct f)
{
  return new real(real x){return y*f(x);};
}

zfct2 operator +(zfct2 f, zfct2 g)
{// Défini f+g
  return new pair(real t, real p){return f(t,p)+g(t,p);};
}

zfct2 operator -(zfct2 f, zfct2 g)
{// Défini f-g
  return new pair(real t, real p){return f(t,p)-g(t,p);};
}

zfct2 operator /(zfct2 f, real x)
{// Défini f/x
  return new pair(real t, real p){return f(t,p)/x;};
}

zfct2 operator *(real x,zfct2 f)
{// Défini x*f
  return new pair(real t, real p){return x*f(t,p);};
}

fct fct(real[] p)
{ // convertit le tableau des coefs du poly p en fonction polynôme
  return new real(real x){
    real y=0;
    for (int i : p.keys) {
      y += p[i]*x^i;
    }
    return y;
  };
}

real C(int l, int m)
{
  if(m < 0) return 1/C(l,-m);
  real OC=1;
  int d=l-m, s=l+m;
  for (int i=d+1; i <=s ; ++i) OC *= i;
  return 1/OC;
}

int csphase=-1;
fct P(int l, int m)
{ // Polynôme de Legendre associé
  // http://mathworld.wolfram.com/LegendrePolynomial.html
  if(m < 0) return (-1)^(-m)*C(l,-m)*P(l,-m);
  real[] xl2;
  for (int k=0; k <= l; ++k) {
    xl2.push((-1)^(l-k)*binomial(l,k));
    if(k != l) xl2.push(0);
  }
  fct dxl2=fct(pdiff(xl2,l+m));
  return new real(real x){
    return (csphase)^m/(2^l*factorial(l))*(1-x^2)^(m/2)*dxl2(x);
  };
}

zfct2 Y(int l, int m)
{// http://fr.wikipedia.org/wiki/Harmonique_sph%C3%A9rique#Expression_des_harmoniques_sph.C3.A9riques_normalis.C3.A9es
  return new pair(real theta, real phi) {
    return sqrt((2*l+1)*C(l,m)/(4*pi))*P(l,m)(cos(theta))*expi(m*phi);
  };
}

real xyabs(triple z){return abs(xypart(z));}

size(16cm);
currentprojection=orthographic(0,1,1);

zfct2 Ylm;

triple F(pair z)
{
  //   real r=0.75+dot(0.25*I,Ylm(z.x,z.y));
  //   return r*expi(z.x,z.y);
  real r=abs(Ylm(z.x,z.y))^2;
  return r*expi(z.x,z.y);
}

int nb=4;
for (int l=0; l < nb; ++l) {
  for (int m=0; m <= l; ++m) {
    Ylm=Y(l,m);

    surface s=surface(F,(0,0),(pi,2pi),60);
    s.colors(palette(s.map(xyabs),Rainbow()));

    triple v=(-m,0,-l);
    draw(shift(v)*s);
    label("$Y_"+ string(l) + "^" + string(m) + "$:",shift(X/3)*v);
  }
}

Étiquettes : , , , , ,


Unofficial package hull_pi.asy – fig0030

Category: Asymptote,hull_pi.asy,Unofficial packagesPh. Ivaldi @ 20 h 09 min

Figure 0003
(Compiled with Asymptote version 2.14svn-r5318)
    
import hull_pi;
import stats;
size(10cm);

pair[] cloud;
int nbpt=200;

// Generate random points.
for (int i=0; i < nbpt; ++i)
  cloud.push((10*unitrand(),10*unitrand()));

////////// GREY PART: depthMax=1 //////////
pair[] hull=hull(cloud,depthMin=0,depthMax=1,angleMin=0,angleMax=360);
filldraw(polygon(hull),lightgrey);

////////// YELLOW PART: depthMax=2 //////////
pair[] hull=hull(cloud,depthMin=0,depthMax=2,angleMin=0,angleMax=360);
filldraw(polygon(hull),lightyellow);

////////// BLUE PART: depthMax=3 //////////
pair[] hull=hull(cloud,depthMin=0,depthMax=3,angleMin=0,angleMax=360);
filldraw(polygon(hull),lightblue);

dot(cloud,red);

Étiquettes :


Unofficial package hull_pi.asy – fig0040

Category: Asymptote,hull_pi.asy,Unofficial packagesPh. Ivaldi @ 21 h 09 min

Figure 0004
(Compiled with Asymptote version 2.14svn-r5318)
    
import hull_pi;
import stats;
size(10cm);

pair[] cloud;
int nbpt=15;
int depthMax=3;

// Generate random points.
for (int i=0; i < nbpt; ++i)
  cloud.push((10*unitrand(),10*unitrand()));

////////// WHITE PART: pivot automaticly computed //////////
pair[] hull=hull(cloud,depthMin=0,depthMax=depthMax,angleMin=0,angleMax=360);
draw(polygon(hull),2mm+white);

// ////////// RED PART: pivot=1 //////////
pair[] hull=hull(cloud,depthMin=0,depthMax=depthMax,angleMin=0,angleMax=360,1);
draw(polygon(hull),1mm+red);

// ////////// BLUE PART: pivot=3 //////////
pair[] hull=hull(cloud,depthMin=0,depthMax=depthMax,angleMin=0,angleMax=360,3);
draw(polygon(hull),bp+blue);

dot(cloud,yellow);
shipout(bbox(3mm,Fill));

Étiquettes :


Official Asymptote example – cheese

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

Figure 0024
(Compiled with Asymptote version 2.14svn-r5318)
/* This code comes from The Official Asymptote Gallery */
    
import graph3;
import palette;
import contour3;
size(400);

real f(real x, real y, real z) {
  return cos(x)*sin(y)+cos(y)*sin(z)+cos(z)*sin(x);
}

surface sf=surface(contour3(f,(-2pi,-2pi,-2pi),(2pi,2pi,2pi),12));
sf.colors(palette(sf.map(abs),Gradient(red,yellow)));
draw(sf,nolight,render(merge=true));

Étiquettes : , , , , ,


Official Asymptote example – condor

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

Figure 0031
(Compiled with Asymptote version 2.14svn-r5318)
/* This code comes from The Official Asymptote Gallery */
    
// Peter Luschny's Condor function
// http://www.luschny.de/math/asy/ElCondorYElGamma.html

import palette;
import graph3;

size(300,300,IgnoreAspect);
currentprojection=orthographic(0,-1,0,center=true);
currentlight=White;
real K=7;

triple condor(pair t)
{
  real y=t.y;
  real x=t.x*y;
  real e=gamma(y+1);
  real ymx=y-x;
  real ypx=y+x;
  real a=gamma((ymx+1)/2);
  real b=gamma((ymx+2)/2);
  real c=gamma((ypx+1)/2);
  real d=gamma((ypx+2)/2);
  real A=cos(pi*ymx);
  real B=cos(pi*ypx);
  return (x,y,log(e)+log(a)*((A-1)/2)+log(b)*((-A-1)/2)+log(c)*((B-1)/2)+
          log(d)*((-B-1)/2));
}

surface s=surface(condor,(-1,0),(1,K),16,Spline);
s.colors(palette(s.map(zpart),Rainbow()));

draw(s,render(compression=Low,merge=true));

É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 – fillcontour

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

Figure 0066
(Compiled with Asymptote version 2.14svn-r5318)
/* This code comes from The Official Asymptote Gallery */
    
import graph;
import palette;
import contour;

size(12cm,IgnoreAspect);

pair a=(pi/2,0);
pair b=(3pi/2,2pi);

real f(real x, real y) {return cos(x)*sin(y);}

int N=100;
int Divs=10;

defaultpen(1bp);

bounds range=bounds(-1,1);
    
real[] Cvals=uniform(range.min,range.max,Divs);
guide[][] g=contour(f,a,b,Cvals,N,operator --);

pen[] Palette=quantize(Rainbow(),Divs);

pen[][] interior=interior(g,extend(Palette,grey,black));
fill(g,interior);
draw(g);

palette("$f(x,y)$",range,point(SE)+(0.5,0),point(NE)+(1,0),Right,Palette,
        PaletteTicks("$%+#0.1f$",N=Divs));

Étiquettes : , , , , ,


Official Asymptote example – gamma3

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

Figure 0076
(Compiled with Asymptote version 2.14svn-r5318)
/* This code comes from The Official Asymptote Gallery */
    
import graph3;
import palette;

size(12cm,IgnoreAspect);
currentprojection=orthographic(1,-2,1);

real X=4.5;
real M=abs(gamma((X,0)));

pair Gamma(pair z) 
{
  return (z.x > 0 || z != floor(z.x)) ? gamma(z) : M;
}

real f(pair z) {return min(abs(Gamma(z)),M);}

surface s=surface(f,(-2.1,-2),(X,2),70,Spline);

real Arg(triple v)
{
  return degrees(Gamma((v.x,v.y)),warn=false);
}

s.colors(palette(s.map(Arg),Wheel()));
draw(s,render(compression=Low,merge=true));

real xmin=point((-1,-1,-1)).x;
real xmax=point((1,1,1)).x;
draw((xmin,0,0)--(xmax,0,0),dashed);

xaxis3("$\mathop{\rm Re} z$",Bounds,InTicks);
yaxis3("$\mathop{\rm Im} z$",Bounds,InTicks(beginlabel=false));
zaxis3("$|\Gamma(z)|$",Bounds,InTicks);

Étiquettes : , , , ,


Official Asymptote example – imagehistogram

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

Figure 0098
(Compiled with Asymptote version 2.14svn-r5318)
/* This code comes from The Official Asymptote Gallery */
    
import stats;
import graph; 
import palette; 
import contour; 
 
size(20cm);

scale(false);

pair[] data=new pair[50000];
for(int i=0; i < data.length; ++i)
  data[i]=Gaussrandpair(); 
 
// Histogram limits and number of bins 
pair datamin=(-0.15,-0.15); 
pair datamax=(0.15,0.15); 
int Nx=30; 
int Ny=30; 

int[][] bins=frequency(data,datamin,datamax,Nx,Ny);
 
real[] values=new real[Nx*Ny]; 
pair[] points=new pair[Nx*Ny];
int k=0; 
real dx=(datamax.x-datamin.x)/Nx;
real dy=(datamax.y-datamin.y)/Ny;
for(int i=0; i < Nx; ++i) {
  for(int j=0; j < Ny; ++j) {
    values[k]=bins[i][j]; 
    points[k]=(datamin.x+(i+0.5)*dx,datamin.y+(j+0.5)*dy); 
    ++k; 
  }
} 
 
// Create a color palette 
pen[] InvGrayscale(int NColors=256) {
  real ninv=1.0/(NColors-1.0); 
  return sequence(new pen(int i) {return gray(1-17*i*ninv);},NColors); 
} 
 
// Draw the histogram, with axes 
bounds range=image(points,values,Range(0,40),InvGrayscale()); 
draw(contour(points,values,new real[] {1,2,3,4,8,12,16,20,24,28,32,36,40},
             operator--),blue); 
xaxis("$x$",BottomTop,LeftTicks,above=true); 
yaxis("$y$",LeftRight,RightTicks,above=true); 


Étiquettes : , , , ,


Official Asymptote example – odetest

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

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

write("integration test");
real f(real t, real x) {return cos(x);}
write(integrate(1,f,0,10,0.1,dynamic=true,0.0002,0.0004,RK3BS,verbose=true));
write();

write("system integration test");
real[] f(real t, real[] x) {return new real[] {x[1],1.5*x[0]^2};}
write(integrate(new real[] {4,-8},f,0,1,n=100,dynamic=true,tolmin=0.0002,
                tolmax=0.0004,RK3BS,verbose=true));
write();

write("simultaneous newton test");
real[] function(real[] x) {
  return new real[] {x[0]^2+x[1]^2-25,(x[0]-6)^2+x[1]^2-25};
}
real[][] fJac(real[] x) {
  return new real[][] {{2*x[0],2*x[1]},{2*(x[0]-6),2*x[1]}};
}
write(newton(function,fJac,new real[] {0,-1}));
write();


write("BVP solver test");
write("Finding initial conditions that solve w''(t)=1.5*w(t), w(0)=4, w(1)=1");
real[] initial(real[] x) {
  return new real[] {4,x[0]};
}

real[] discrepancy(real[] x) {
  write("Error: ",x[0]-1);
  return new real[] {x[0]-1};
}

write(solveBVP(f,0,1,n=10,initial,discrepancy,guess=new real[] {-30},RK4,
               iterations=10));
write();
write(solveBVP(f,0,1,n=100,initial,discrepancy,guess=new real[] {-30},RK4,
               iterations=10));
write();
write(solveBVP(f,0,1,n=10000,initial,discrepancy,guess=new real[] {-30},RK4,
               iterations=10));
write();

Étiquettes :


Official Asymptote example – spectrum

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

Figure 0209
(Compiled with Asymptote version 2.14svn-r5318)
/* This code comes from The Official Asymptote Gallery */
    
import graph;
usepackage("ocg");
settings.tex="pdflatex";

// Dan Bruton algorithm
pen nm2rgb(real wl, real gamma=0.8, bool intensity=true) {
  triple rgb;
  if(wl >= 380 && wl <= 440) {rgb=((440-wl)/60,0,1);}
  if(wl >  440 && wl <= 490) {rgb=(0,(wl-440)/50,1);}
  if(wl >  490 && wl <= 510) {rgb=(0,1,(510-wl)/20);}
  if(wl >  510 && wl <= 580) {rgb=((wl-510)/70,1,0);}
  if(wl >  580 && wl <= 645) {rgb=(1,(645-wl)/65,0);}
  if(wl >  645 && wl <= 780) {rgb=(1,0,0);}
  
  real Intensity=1;
  if(intensity) {
    if(wl >= 700) {Intensity=0.3+0.7*(780-wl)/80;}
    else if(wl <= 420) {Intensity=0.3+0.7*(wl-380)/40;}
  }

  return rgb((Intensity*rgb.x)**gamma,(Intensity*rgb.y)**gamma,
             (Intensity*rgb.z)**gamma);
}

real width=1;
real height=50;

begin("spectrum");
for(real i=380 ; i <= 780 ; i += width) {
  draw((i,0)--(i,height),width+nm2rgb(wl=i,false)+squarecap);
}
begin("Extinction",false); // nested
for(real i=380 ; i <= 780 ; i += width) {
  draw((i,0)--(i,height),width+nm2rgb(wl=i,true)+squarecap);
}
end();
end();

begin("Wavelength");
xaxis(scale(0.5)*"$\lambda$(nm)",BottomTop,380,780,
      RightTicks(scale(0.5)*rotate(90)*Label(),step=2,Step=10),above=true);
end();

// From Astronomical Data Center(NASA)
// Neutral only
real[] Na={423.899, 424.208, 427.364, 427.679, 428.784, 429.101,
           432.14, 432.462, 434.149, 434.474, 439.003, 439.334, 441.989, 442.325,
           449.418, 449.766, 454.163, 454.519, 568.2633, 568.8204, 588.995,
           589.5924};
begin("Na absorption");
for(int i=0; i < Na.length; ++i) {
  draw((Na[i],0)--(Na[i],height),0.1*width+squarecap);
}
end();

begin("Na emission");
for(int i=0; i < Na.length; ++i) {
  draw((Na[i],0)--(Na[i],-height),0.1*width+nm2rgb(Na[i],false)+squarecap);
}
end();

// Neutral only
real[] Zn={388.334, 396.543, 411.321, 429.288, 429.833, 462.981,
           468.014, 472.215, 481.053 , 506.866, 506.958, 518.198, 530.865,
           531.024, 531.102, 577.21, 577.55, 577.711, 623.79, 623.917, 636.234,
           647.918, 692.832, 693.847, 694.32, 779.936};
begin("Zn absorption",false);
for(int i=0; i < Zn.length; ++i) {
  draw((Zn[i],0)--(Zn[i],height),width+squarecap);
}
end();

begin("Zn emission",false);
for(int i=0; i < Zn.length; ++i) {
  draw((Zn[i],0)--(Zn[i],-height),width+nm2rgb(Zn[i],false)+squarecap);
}
end();

shipout(bbox(2mm,Fill(white)));

Étiquettes : , , , ,


Official Asymptote example – triangulate

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

Figure 0241
(Compiled with Asymptote version 2.14svn-r5318)
/* This code comes from The Official Asymptote Gallery */
    
size(200);
int np=100;
pair[] points;

real r() {return 1.2*(rand()/randMax*2-1);}

for(int i=0; i < np; ++i)
  points.push((r(),r()));

int[][] trn=triangulate(points);

for(int i=0; i < trn.length; ++i) {
  draw(points[trn[i][0]]--points[trn[i][1]]);
  draw(points[trn[i][1]]--points[trn[i][2]]);
  draw(points[trn[i][2]]--points[trn[i][0]]);
}

for(int i=0; i < np; ++i)
  dot(points[i],red);

Étiquettes : ,


Animation with Asymptote – fig0050

Category: Animation,AsymptotePh. Ivaldi @ 12 h 20 min

Figure 0005
(Compiled with Asymptote version 1.43)
Movie flash (swf)
This animation is available in the Syracuse web site.
    
/* Author: Nathan Carter. */
include "./makecd.asy";
import animate;
// settings.tex="pdflatex";
settings.keep=true;

animation A;
A.global=false;

real length = 4; // seconds
int fps = 50;
real rad = 6;
real ht = 2;
real pixsz = 300;
real ptsz = 1;
int loops = 6;
pen border = black;

real frames = length*fps;
picture tmp;

size(pixsz);
for (int i=100 ; i < 100+frames ; ++i) {
  save();
  add(CayleyDiagram(nodeLocs, arrows, orders, arrowPens,
                    cam = (rad*cos(2*i*pi/frames),rad*sin(2*i*pi/frames),ht),
                    arrowThickness = 2, nodeSize = 0.02,
                    arrowMargin = 1mm, depthCueing = true ));
  draw(box((-ptsz/2,-ptsz/2), (ptsz/2,ptsz/2)), border);
  A.add();
  write( "Did " + (string)(i-99) + " out of " + (string)frames );
  restore();
}

write( "Merging..." );
A.movie(delay=(int)(100/fps));
write( "Done." );

Étiquettes : ,


Animation with Asymptote – fig0190

Category: Animation,AsymptotePh. Ivaldi @ 2 h 20 min

Figure 0019
(Compiled with Asymptote version 1.86svn-r4626)
Movie flash (swf)
This animation is available in the Syracuse web site.
    
size(12cm,0);
import hull_pi;
import animation;
animation A;
settings.outformat="pdf";

pair[] cloud;
int nbpt=200;

// Generate random points.
for (int i=0; i < nbpt; ++i)
  cloud.push((10*unitrand(),10*unitrand()));

for (int i=1; i < 8; ++i) { // Control the depth
  for (int j=0; j < 30; ++j) { // Point of view = cloud[i]
    real depthMax=i/2;
    // Nodes of the hull
    pair[] hull=hull(cloud,depthMin=0,depthMax=depthMax,angleMin=10,angleMax=360,pivot=j);

    save();// Add new picture to the animation
    filldraw(polygon(hull),lightgrey);
    dot(cloud[j],3mm+green);
    dot(cloud,red);
    label("depthMax="+(string)depthMax);
    A.add();
    restore();
  }
}

A.movie();

Étiquettes : ,