メモ: わからないながらにPlaskで遊んでみた2

結果:

Plask:

http://www.plask.org/

参考:

web雑記帳: javascriptで4次Runge-Kutta法

重振り子

コード:

var plask = require('plask');

//
// refer to:
// http://yoshifuji.blogspot.jp/2011/07/javascript4runge-kutta.html
// http://www.ne.jp/asahi/tokyo/nkgw/gakusyu/rikigaku/2zyufuriko/zyufuriko_kaisetu/zyufuriko_kaisetu.html
function solver(x1,x2,x3,x4,t,dt){  
  var k11,k12,k13,k14,k21,k22,k23,k24,k31,k32,k33,k34,k41,k42,k43,k44;  
  
  k11 = dt * func1(x1,x2,x3,x4,t);  
  k12 = dt * func2(x1,x2,x3,x4,t);  
  k13 = dt * func3(x1,x2,x3,x4,t);  
  k14 = dt * func4(x1,x2,x3,x4,t);  
   
  k21 = dt * func1(x1 + k11/2, x2 + k12/2, x3 + k13/2, x4 + k14/2, t + dt/2);  
  k22 = dt * func2(x1 + k11/2, x2 + k12/2, x3 + k13/2, x4 + k14/2, t + dt/2);  
  k23 = dt * func3(x1 + k11/2, x2 + k12/2, x3 + k13/2, x4 + k14/2, t + dt/2);  
  k24 = dt * func4(x1 + k11/2, x2 + k12/2, x3 + k13/2, x4 + k14/2, t + dt/2);  

  k31 = dt * func1(x1 + k21/2, x2 + k22/2, x3 + k23/2, x4 + k24/2, t + dt/2);  
  k32 = dt * func2(x1 + k21/2, x2 + k22/2, x3 + k23/2, x4 + k24/2, t + dt/2);  
  k33 = dt * func3(x1 + k21/2, x2 + k22/2, x3 + k23/2, x4 + k24/2, t + dt/2);  
  k34 = dt * func4(x1 + k21/2, x2 + k22/2, x3 + k23/2, x4 + k24/2, t + dt/2);  

  k41 = dt * func1(x1 + k31, x2 + k32, x3 + k33/2, x4 + k34/2, t + dt);  
  k42 = dt * func2(x1 + k31, x2 + k32, x3 + k33/2, x4 + k34/2, t + dt);  
  k43 = dt * func3(x1 + k31, x2 + k32, x3 + k33/2, x4 + k34/2, t + dt);  
  k44 = dt * func4(x1 + k31, x2 + k32, x3 + k33/2, x4 + k34/2, t + dt);  

  var tmp= [];
  tmp[0] = x1 + (k11/6 + k21/3 + k31/3 + k41/6);  
  tmp[1] = x2 + (k12/6 + k22/3 + k32/3 + k42/6);  
  tmp[2] = x3 + (k13/6 + k23/3 + k33/3 + k43/6);  
  tmp[3] = x4 + (k14/6 + k24/3 + k34/3 + k44/6);

  return tmp;  
}
var m1= 1.;
var m2= 2.;
var l1= 9.8;
var l2= 9.8;
var  g= 9.8; //...
function func1(x1,x2,x3,x4,t){  
  return x2*0.99999;
}  
function func2(x1,x2,x3,x4,t){
  var c0= x1-x3;
  var cc= Math.cos(c0);
  var cs= Math.sin(c0);

  var c1= 1.0/(m1+m2*cs*cs);
  var c2= -(m1+m2)*Math.sin(x1)/l1*g;
  var c3=  m2*cc*Math.sin(x2)/l1*g;
  var c4= -m2*cc*cs*x2*x2;
  var c5= -m2*cs*(l2/l1)*x4*x4;
  return c1*(c2+c3+c4+c5); 
}  
function func3(x1,x2,x3,x4,t){  
  return x4*0.99999;
}
function func4(x1,x2,x3,x4,t){  
  var c0= x1-x3;
  var cc= Math.cos(c0);
  var cs= Math.sin(c0);
  
  var c1= 1.0/(m1+m2*cs*cs);
  var c2= (m1+m2)*cs*Math.cos(x1)/l2*g;
  var c3= (m1+m2)*cs*Math.sin(x2)*(l1/l2)*x2*x2;
  var c4= m2*cc*cs*x4*x4;
  return c1*(c2+c3+c4);
}  
//

var path= new PlaskRawMac.SkPath;

var dots1= [];
var dots2= [];
var dotsN= 90.;

var xs = [];
var t;
var dt;

var flag;
plask.simpleWindow({
  settings: {
    width: 400,
    height: 400
  },

  init: function() {
    var canvas = this.canvas, paint = this.paint;
    this.framerate(60);

    flag= false;
    this.on('keyDown', function(e) {
      switch (e.str) {
        case ' ': flag= true; break;
      }
    });
    this.on('keyUp', function(e) {
      switch (e.str) {
        case ' ': flag= false; break;
      }
    });

    for(var i=0; i<dotsN; i++) {
      dots1.push({x:0, y:0});
      dots2.push({x:0, y:0});
    }
    xs[0]= Math.PI/3;
    xs[1]= 0.1;
    xs[2]= Math.PI;
    xs[3]= -0.1;
    dt= 1./720.;
    t= 0.0;
  },

  draw: function() {
    var canvas = this.canvas, paint = this.paint;
    canvas.clear(0, 0, 0, 255);    

    xs= solver(xs[0], xs[1], xs[2], xs[3], t, dt);
    t+= dt;
    xs= solver(xs[0], xs[1], xs[2], xs[3], t, dt);
    t+= dt;
    xs= solver(xs[0], xs[1], xs[2], xs[3], t, dt);
    t+= dt;
    xs= solver(xs[0], xs[1], xs[2], xs[3], t, dt);
    t+= dt;

    var y= 200;
    var x= 200;
    var th1= xs[0];
    var th2= xs[2];
    y+= 10.*l1*Math.cos(th1);
    x+= 10.*l1*Math.sin(th1);
    dots1.push({x:x, y:y});
    dots1.shift();
    y+= 10.*l2*Math.cos(th2);
    x+= 10.*l2*Math.sin(th2);
    dots2.push({x:x, y:y});
    dots2.shift();

    xs= solver(xs[0], xs[1], xs[2], xs[3], t, dt);
    t+= dt;
    xs= solver(xs[0], xs[1], xs[2], xs[3], t, dt);
    t+= dt;
    xs= solver(xs[0], xs[1], xs[2], xs[3], t, dt);
    t+= dt;
    xs= solver(xs[0], xs[1], xs[2], xs[3], t, dt);
    t+= dt;

    y= 200;
    x= 200;
    var th1= xs[0];
    var th2= xs[2];
    y+= 10.*l1*Math.cos(th1);
    x+= 10.*l1*Math.sin(th1);
    dots1.push({x:x, y:y});
    dots1.shift();
    y+= 10.*l2*Math.cos(th2);
    x+= 10.*l2*Math.sin(th2);
    dots2.push({x:x, y:y});
    dots2.shift();


    xs= solver(xs[0], xs[1], xs[2], xs[3], t, dt);
    t+= dt;
    xs= solver(xs[0], xs[1], xs[2], xs[3], t, dt);
    t+= dt;
    xs= solver(xs[0], xs[1], xs[2], xs[3], t, dt);
    t+= dt;
    xs= solver(xs[0], xs[1], xs[2], xs[3], t, dt);
    t+= dt;

    path.reset();
    y= 200;
    x= 200;
    path.moveTo(x, y);
    var th1= xs[0];
    var th2= xs[2];
    y+= 10.*l1*Math.cos(th1);
    x+= 10.*l1*Math.sin(th1);
    path.lineTo(x, y);
    dots1.push({x:x, y:y});
    dots1.shift();
    y+= 10.*l2*Math.cos(th2);
    x+= 10.*l2*Math.sin(th2);
    path.lineTo(x, y);
    dots2.push({x:x, y:y});
    dots2.shift();


    if (flag) {
      paint.setStroke();
      paint.setColor(125,125,125,200);
      canvas.drawPath(paint, path);
    }

    paint.setFill();
    for(var i= 0; i<dotsN; i++) {
      var p1= dots1[i];
      paint.setColor(255, 125, 125, 255/dotsN*i);
      canvas.drawCircle(paint, p1.x, p1.y, i/dotsN*10.);
      var p2= dots2[i];
      paint.setColor(125, 125, 355, 255/dotsN*i);
      canvas.drawCircle(paint, p2.x, p2.y, i/dotsN*10.);
    }
  }
});

メモ: わからないながらにPlaskで遊んでみた

結果:

Plask:

http://www.plask.org/

参考:

web雑記帳: javascriptで4次Runge-Kutta法

重振り子

コード:

var plask = require('plask');

//
// refer to:
// http://yoshifuji.blogspot.jp/2011/07/javascript4runge-kutta.html
// http://www.ne.jp/asahi/tokyo/nkgw/gakusyu/rikigaku/2zyufuriko/zyufuriko_kaisetu/zyufuriko_kaisetu.html
function solver(x1,x2,x3,x4,t,dt){  
  var k11,k12,k13,k14,k21,k22,k23,k24,k31,k32,k33,k34,k41,k42,k43,k44;  
  
  k11 = dt * func1(x1,x2,x3,x4,t);  
  k12 = dt * func2(x1,x2,x3,x4,t);  
  k13 = dt * func3(x1,x2,x3,x4,t);  
  k14 = dt * func4(x1,x2,x3,x4,t);  
   
  k21 = dt * func1(x1 + k11/2, x2 + k12/2, x3 + k13/2, x4 + k14/2, t + dt/2);  
  k22 = dt * func2(x1 + k11/2, x2 + k12/2, x3 + k13/2, x4 + k14/2, t + dt/2);  
  k23 = dt * func3(x1 + k11/2, x2 + k12/2, x3 + k13/2, x4 + k14/2, t + dt/2);  
  k24 = dt * func4(x1 + k11/2, x2 + k12/2, x3 + k13/2, x4 + k14/2, t + dt/2);  

  k31 = dt * func1(x1 + k21/2, x2 + k22/2, x3 + k23/2, x4 + k24/2, t + dt/2);  
  k32 = dt * func2(x1 + k21/2, x2 + k22/2, x3 + k23/2, x4 + k24/2, t + dt/2);  
  k33 = dt * func3(x1 + k21/2, x2 + k22/2, x3 + k23/2, x4 + k24/2, t + dt/2);  
  k34 = dt * func4(x1 + k21/2, x2 + k22/2, x3 + k23/2, x4 + k24/2, t + dt/2);  

  k41 = dt * func1(x1 + k31, x2 + k32, x3 + k33/2, x4 + k34/2, t + dt);  
  k42 = dt * func2(x1 + k31, x2 + k32, x3 + k33/2, x4 + k34/2, t + dt);  
  k43 = dt * func3(x1 + k31, x2 + k32, x3 + k33/2, x4 + k34/2, t + dt);  
  k44 = dt * func4(x1 + k31, x2 + k32, x3 + k33/2, x4 + k34/2, t + dt);  

  var tmp= [];
  tmp[0] = x1 + (k11/6 + k21/3 + k31/3 + k41/6);  
  tmp[1] = x2 + (k12/6 + k22/3 + k32/3 + k42/6);  
  tmp[2] = x3 + (k13/6 + k23/3 + k33/3 + k43/6);  
  tmp[3] = x4 + (k14/6 + k24/3 + k34/3 + k44/6);

  return tmp;  
}  
function func1(x1,x2,x3,x4,t){  
  return x2;  
}  
function func2(x1,x2,x3,x4,t){  
  return -2 * x1 + 1 * x3; 
}  
function func3(x1,x2,x3,x4,t){  
  return x4;
}
function func4(x1,x2,x3,x4,t){  
  return 2 * x1 - 2 * x3;
}  
//

var path= new PlaskRawMac.SkPath;

var dots1= [];
var dots2= [];
var dotsN= 30.;

var xs = [];
var t;
var dt;

var flag;
plask.simpleWindow({
  settings: {
    width: 400,
    height: 400
  },

  init: function() {
    var canvas = this.canvas, paint = this.paint;
    this.framerate(60);

    flag= false;
    this.on('keyDown', function(e) {
      switch (e.str) {
        case ' ': flag= true; break;
      }
    });
    this.on('keyUp', function(e) {
      switch (e.str) {
        case ' ': flag= false; break;
      }
    });

    for(var i=0; i<dotsN; i++) {
      dots1.push({x:0, y:0});
      dots2.push({x:0, y:0});
    }
    xs[0]= Math.PI;
    xs[1]= 1.0;
    xs[2]= Math.PI;
    xs[3]=-5.0;
    dt= 1./180.;
    t= 0.0;
  },

  draw: function() {
    var canvas = this.canvas, paint = this.paint;
    canvas.clear(0, 0, 0, 255);    

    xs= solver(xs[0], xs[1], xs[2], xs[3], t, dt);
    t+= dt;
    xs= solver(xs[0], xs[1], xs[2], xs[3], t, dt);
    t+= dt;
    xs= solver(xs[0], xs[1], xs[2], xs[3], t, dt);
    t+= dt;

    path.reset();
    var y= 200;
    var x= 200;
    path.moveTo(x, y);
    var th1= xs[0];
    var th2= xs[2];
    y+= 98*Math.cos(th1);
    x+= 98*Math.sin(th1);
    path.lineTo(x, y);
    dots1.push({x:x, y:y});
    dots1.shift();
    y+= 98*Math.cos(th2);
    x+= 98*Math.sin(th2);
    path.lineTo(x, y);
    dots2.push({x:x, y:y});
    dots2.shift();


    if (flag) {
      paint.setStroke();
      paint.setColor(125,125,125,200);
      canvas.drawPath(paint, path);
    }

    paint.setFill();
    for(var i= 0; i<dotsN; i++) {
      var p1= dots1[i];
      paint.setColor(255, 125, 125, 255/dotsN*i);
      canvas.drawCircle(paint, p1.x, p1.y, i/dotsN*10.);
      var p2= dots2[i];
      paint.setColor(125, 125, 355, 255/dotsN*i);
      canvas.drawCircle(paint, p2.x, p2.y, i/dotsN*10.);
    }
  }
});

Quartz Composer / GLSL Shader / BokehなFragment Shaderを生成するコード

参考

THE ART-LOG OF MARTINS UPITIS: a GLSL depth of field filter with bokeh

GLSL depth of field with bokeh v2.4 (update)

コード
参考のコード他をよくわからないまま弄くり回していました．
とりあえずサンプリングする点を求めるのにコードにしたのでそれをメモしておきます．

#include <stdio.h>
#include <stdlib.h>
double *cross(double *p, double *q) {
 double *r;
 r= (double *)malloc(3*sizeof(double));
 r[0]= p[1]*q[2]-p[2]*q[1];
 r[1]=-p[0]*q[2]+p[2]*q[0];
 r[2]= p[0]*q[1]-p[1]*q[0];
 return r;
}
double *sub(double *p, double *q) {
 double *r;
 r= (double *)malloc(3*sizeof(double));
 r[0]= p[0]-q[0];
 r[1]= p[1]-q[1];
 r[2]= p[2]-q[2];
 return r;
}
int isinCircle(double *q) {
 return (q[0]*q[0]+q[1]*q[1]<0.25)?1:0;
}
int isin6(double *q)
{
 double p[6][3]= {
     { 0.0, 0.50, 0.},
     { 0.5, 0.25, 0.},
     { 0.5,-0.25, 0.},
     { 0.0,-0.50, 0.},
     {-0.5,-0.25, 0.},
     {-0.5, 0.25, 0.}};
 double *r;
 int i, c= 0;

 for(i= 0; i<5; i++) {
  r= cross(sub(p[i], q), sub(p[i+1], q)); 
  c+= r[2]<0?0:1;
 }
 r= cross(sub(p[5], q), sub(p[0], q)); 
 c+= r[2]<0?0:1;

 return (c==0)?1:0;
}
int main(int argc, char* argv[]) {

 double q[3]= {0., 0., 0.};
 int j,k,c= 0;
 int r= 5;
 
 printf("uniform sampler2D RenderedTexture;\nuniform sampler2D DepthTexture;\n");
 printf("const float blurclamp = 1.;\nconst float bias = 10.;\nuniform float focus;\n");
 printf("uniform float pixelsWide;\nuniform float pixelsHigh;\n");
 printf("vec2 texcel = vec2(1./pixelsWide, 1./pixelsHigh);\n");
 printf("void main()\n{\n");
 printf("\tvec4 depth = texture2D(DepthTexture,gl_TexCoord[0].xy );\n");
 printf("\tfloat factor = ( depth.x - focus );\n");
 printf("\tvec2 dofblur = vec2 (clamp( factor * bias, -blurclamp, blurclamp ));\n");
 printf("\tvec4 col = vec4(0.0);\n");
 
 for(k= -r; k<r; k++) {
  for(j= -r; j<r; j++) {
   q[0]= k/(r*2.);
   q[1]= j/(r*2.);
   if (isin6(q)==1) {
    if (0) {
     printf("\tcol += texture2D(RenderedTexture, gl_TexCoord[0].xy + (vec2(");
     printf("%2d., %2d.", k, j);
     printf(")*texcel) * dofblur);\n");
    }
    else {
     printf("\tcol.r += texture2D(RenderedTexture, gl_TexCoord[0].xy + (vec2(");
     printf("%2d., %2d.", k, j);
     printf(")*texcel) * dofblur*vec2( .000, .50)).r;\n");
     printf("\tcol.g += texture2D(RenderedTexture, gl_TexCoord[0].xy + (vec2(");
     printf("%2d., %2d.", k, j);
     printf(")*texcel) * dofblur*vec2( .866,-.25)).g;\n");
     printf("\tcol.b += texture2D(RenderedTexture, gl_TexCoord[0].xy + (vec2(");
     printf("%2d., %2d.", k, j);
     printf(")*texcel) * dofblur*vec2(-.866,-.25)).b;\n");
    }
    c++;
   }
  }
 }
 
 printf("\tgl_FragColor = col/%d.;\n\tgl_FragColor.a = 1.0;\n}\n", c);
 return 0;
}

結果
こんな感じ．