Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] Diff of /branches/pure-cfg/test/lic.diderot
ViewVC logotype

Diff of /branches/pure-cfg/test/lic.diderot

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 935, Fri Apr 22 13:22:41 2011 UTC revision 937, Fri Apr 22 13:49:34 2011 UTC
# Line 5  Line 5 
5  //  //
6  // process output with:  // process output with:
7  // unu reshape -i mip.txt -s 300 200 | unu quantize -b 8  -o lic.png  // unu reshape -i mip.txt -s 300 200 | unu quantize -b 8  -o lic.png
 //  
 // NOTE: this computes a streamline along only one direction of  
 // field at starting point; it would be more correct to go  
 // in both directions (both upstream and downstream)  
8    
9  int imgSizeX = 300;  int imgSizeX = 300;
10  int imgSizeY = 200;  int imgSizeY = 200;
11  real h = 0.01;  real h = 0.01;
12  int stepNum = 20;  int stepNum = 25;
13    
14  image(2)[2] Vimg = load("../data/vorttest.nrrd");  image(2)[2] Vimg = load("../data/vorttest.nrrd");
15  field#1(2)[2] V = Vimg ⊛ ctmr;  field#1(2)[2] V = Vimg ⊛ ctmr;
# Line 24  Line 20 
20      real xx = lerp(0.0, 3.0, -0.5, real(xi), real(imgSizeX)-0.5);      real xx = lerp(0.0, 3.0, -0.5, real(xi), real(imgSizeX)-0.5);
21      real yy = lerp(0.0, 2.0, -0.5, real(yi), real(imgSizeY)-0.5);      real yy = lerp(0.0, 2.0, -0.5, real(yi), real(imgSizeY)-0.5);
22      vec2 pos0 = [xx,yy];      vec2 pos0 = [xx,yy];
23      vec2 pos = pos0;      vec2 forw = pos0;
24      output real avg = 0.0;      vec2 back = pos0;
25        output real avg = R(pos0);
26      int step = 0;      int step = 0;
27    
28      update {      update {
29          avg = avg + R(pos);          // Euler integration step
30          // pos = pos + h*V(pos);         // Euler step          // forw = forw + h*V(forw);
31          pos = pos + h*V(pos + 0.5*h*V(pos));  // RK2 step          // back = back - h*V(back);
32            // Midpoint method step
33            forw = forw + h*V(forw + 0.5*h*V(forw));
34            back = back - h*V(back - 0.5*h*V(back));
35            avg = avg + R(forw) + R(back);
36          step = step + 1;          step = step + 1;
37          if (step == stepNum) {          if (step == stepNum) {
38            // modulate output by velocity at initial position            // modulate output by velocity at initial position
39            avg = |V(pos0)|*avg/real(stepNum);            avg = |V(pos0)|*avg/real(1 + 2*stepNum);
40            stabilize;            stabilize;
41          }          }
42      }      }

Legend:
Removed from v.935  
changed lines
  Added in v.937

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0