windo = scf(0);         // set current figure (window #0)
windo.color_map = jetcolormap(64);       // color-map new style ???
// xset("colormap",jetcolormap(64));     // color-map old style
// isoview(-1,1,-1,1); // orthogonale diagramme - old style

clf(0);             // clear current figure
drawlater();        // wait for plot complete
clc();              // clear command window
clear();            // kill old variables

 printf("\n\nrefractory waves, 2 dimensional\n");
 printf("example with gaussian waves, euklidic distance\n");
 printf("and homogeneous wave space (not true for neuropile)\n");
 printf("written with Scilab 5.4.1 under WindowsXP\n");
 printf("source www.gheinz.de/publications/animations/index.htm\n");
 mkdir('test'); 
 printf("creates subdirectory \\test for GIF-images in directory \n\n");
 dire=pwd(); datum = date();

// inits /////////////////////////////////////////////

     // physical parameters:
     xmin=-14; xdelta=.5; xmax=14;       // length range x  [for example in mm]
     ymin=-12; ydelta=.5; ymax=12;       // length range y  [for example in mm]
     vtmin=-10; vtdelta=.5; vtmax=30;   // length range vt [v*t = (mm/s)*ms = µm]
     zmin=-1; zmax=4;                   // amplitude range
     refract=10;               // length of refractority period "burned region"
     
     // waves:
     r1=4; x1=10; y1=7;       // wave 1 origin
     r2=2; x2=10; y2=-7;      // wave 2 origin
     r3=9; x3=-10; y3=7;      // wave 3 origin
     r4=7; x4=-10; y4=-7;     // wave 4 origin
     
     // titel
     titelw='waves with valleys';       // waves
     titeli='grassfire waves with refractory phase';                             // integral
     source="http://www.gheinz.de";
     giffiles='test\grassfire4w';
     
     // images and pixels:
     ix = 1;    // x  index (lowest possible is one)
     iy = 1;    // y  index (lowest possible is one)
     ivt = 1;   // vt index (lowest possible is one)
     j=1000;    // gif-image index from 1000, else sorting problems
     box_x=-7; box_y=-11; // vt-label without boxframe

    // function definitions:
 deff('[y] = gauss(u)','y = exp(-(u)^2)');                  // Gaussian wave
 deff('[z] = welle(w)','z = gauss(w) - gauss(w-2)');      // wave with depth
 // deff('y=tal(u)','if u<0 | u>refract; y=0; else y=-zmin; end');  // square function
 // deff('[z] = welle(w)','z = gauss(w) - tal(w-2)');       // wave with depth

// plots ////////////////////////////////////////////////

  xm = xmin:xdelta:xmax; // from:delta:to, vector for x, image size
  ym = ymin:ydelta:ymax; // from:delta:to, vextor for y, image size
      
  for vt = vtmin:vtdelta:vtmax;    // from : delta : to    parameter is vt
          drawlater(); clf(0);     // clear figure window #0
          
          // Alle Pixel des Feldes berechnen:
          for x = xmin:xdelta:xmax;     // from:delta:to
             for y = ymin:ydelta:ymax;  // from:delta:to
               
               // wave values (z):
               welle1 = welle(vt + r1 - sqrt((x-x1)^2 + (y-y1)^2));      // 
               welle2 = welle(vt + r2 - sqrt((x-x2)^2 + (y-y2)^2));      // 
               welle3 = welle(vt + r3 - sqrt((x-x3)^2 + (y-y3)^2));      // 
               welle4 = welle(vt + r4 - sqrt((x-x4)^2 + (y-y4)^2));      // 
               
               // left window, superposition
               ff(ix,iy) = welle1 + welle2 + welle3 + welle4;  // size (xm x ym)
               
               // right window, burning grassfires:
               ii(ix,iy) = min(welle1,welle2,welle3,welle4);
               
               // upper limitation
               if ii(ix,iy)<0 then ii(ix,iy)=zmin; end

               iy = iy+1;
             end;             
             iy = 1;     // loopend, reset
             ix = ix +1;
          end;
          ix = 1;     // loopend, reset
          
          // reset for a new image
          if ivt == 1 then gg = 0; end;
          
          // integration (summation), image matrix size (xm' x ym)          
          gg = min(gg, ii) + ff;     // matrix addition
          
          // add the burned zone to wave field
          // gg = gg + ff;     // matrix addition
          
        // wave plot for actual vt:
          subplot(121);     // Reihen, Spalten, Element
          Sgrayplot(xm,ym,ff,frameflag=4,zminmax=[zmin,zmax]);
          // xset("colormap",jetcolormap(64));  // old style
          //colorbar(zmin,zmax);
          xtitle(titelw,'x','y');
          xstring(box_x,-box_y-1,' vt = '+string(vt)+' ');
          xstring(box_x,box_y,source);             
          
        // integral plot for actual vt:
          subplot(122);     // Reihen, Spalten, Element
          Sgrayplot(xm,ym,gg,frameflag=4,zminmax=[zmin,zmax]);
          // xset("colormap",jetcolormap(64));  // old style
          // colorbar(min(gg),max(gg));
          //colorbar(zmin,zmax);
          xtitle(titeli,'x','y');           // Titel
          xstring(box_x,box_y,source);             

          // xstring(box_x+2,-box_y-1,datum); 
          drawnow();    // jetzt erst zeichnen
          show_window(0);  // window on top, new version
          // xselect();   // window on top, old version
       
        // gif saving:
          gifname = giffiles+string(j)+'.gif';
          xs2gif(0,gifname);                 // save #0 output as *.GIF-file
          printf(" saved:   %s\n",gifname);
          j = j+1;        // index-counter for GIF-images
          ivt = ivt + 1;  // slope-index counter
  end;  // vt; 
printf("\nsaved at \n\n-->"); pwd()   // Konsolenausgabe
return;
