mc.ini

! freopen(stderr,"/dev/null","w");

function clear
{
   "\e[H\e[2J";
!  "\eH\eJ";
} 

!clear
"                         MATC Matrix Language                              "
"                            (Version 0.9)                                  "
" "

function tensor(x,y,oper)
{
   if ((~exists("x")) | (~exists("y")))
   {
     "Usage: tensor(x,y[,oper])";
   } else 
   { 
     if (~exists("oper")) oper = "*";

     make = "tmp = x[i,j]" oper "y";

     nx = size(x)[0]-1;
     mx = size(x)[1]-1;
     ny = size(y)[0];
     my = size(y)[1];
     k  = 0:ny-1;

     for(i=0:nx)
     {
       l  = 0:my-1;
       for(j=0:mx)
       {
         rubbish = @make;
         _tensor[k,l] = tmp;
         l = l + my;
       }
       k = k + ny;
     }
   }
}

function pointw(x,y,oper)
{
   if ((~exists("x")) | (~exists("y")))
   {
       "Usage: pointw(x,y[,oper])";
   } else 
   {
       if (~exists("oper")) oper = "*";

       make = "tmp = (x[i,j])" oper "(y[i,j])";

       nx = size(x)[0]-1;
       mx = size(x)[1]-1;

       for(i=0:nx)
       {
           for(j=0:mx)
           {
               rubbish = @make;
               _pointw[i,j] = tmp;
           }
       }
    }
}

function vdot(x,y)
!
! Return dot product of two arrays of vectors x(0:2,n) and y(0:2,n)
!
{
    n = size(x)[1]
    n = 0:n-1 

    s = x*y
    _vdot(n) = s(0,n) + s(1,n) + s(2,n)
}

function vcross(x,y)                                   
!
! Return cross product of two arrays of vectors x(0:2,n) and y(0:2,n)
!
{
    n = size(x)[1] 
    n = 0:n-1       

    _vcross(2,n) = x(0,n)*y(1,n) - x(1,n)*y(0,n);
    _vcross(1,n) = x(2,n)*y(0,n) - x(0,n)*y(2,n);
    _vcross(0,n) = x(1,n)*y(2,n) - x(2,n)*y(1,n);
}

function time(t,x)
!
! Return indexes for accessing given timestep of variables of current
! element model loaded in elmer-post.
!
! for example: t=velo(0,time(5)) gives x-component of velocity in sixth 
! timestep (counting begins from zero). time("all") gives indexes for
! the whole timestep range.
!
    import temp,nodes
{
    if ( t=="all" )
    {
        if ( exists("x") ) {
          _time = 0:size(x)[1]-1;
        } else {
          "Usage: time(\"all\", x);"
        }
    } else {
        n = size(nodes)[1]
        _time = t*n + (0:(n-1))
    }
}


function dfdt(f,dt)
!
!  compute time derivative of given scalar or vector function
!
   import nodes
{
   n  = size(nodes)[1];
   nf = size(f)[0];
   nt = size(f)[1] / n;

   if ( ~exists("dt") ) dt = 1

   _dfdt = 0;
   ddt = 1.0 / dt

   if ( nt > 1 )
   {
      k0 = time(0);
      k1 = time(1);
      for( j=0:nf-1 )
      {
         _dfdt(j,k0) = ddt * ( f(j,k1) - f(j,k0) )
      }

      for( k=1:nt-2 )
      {
         k0 = time(k-1);
         k1 = time(k);
         k2 = time(k+1);
         for( j=0:nf-1 )
         {
            _dfdt(j,k1) = 0.5 * ddt * ( f(j,k2) - f(j,k0) )
         }
      }

      k0 = time(nt-2);
      k1 = time(nt-1);
      for( j=0:nf-1 )
      {
         _dfdt(j,k1) = ddt * ( f(j,k1) - f(j,k0) )
      }
   }
}

function norm(x,t)
export iter,eps
import eps,Inf
{
   if (~exists("x"))
   {
     "Usage: norm(x[,t])";
   } else 
   {
     if (~exists("eps")) eps = 1.0e-16;
     if (~exists("t")) t = 2;

     rows = size(x)[0];
     cols = size(x)[1];

     if (rows > 1 & cols > 1)
     {
       if (t == 1) 
       {
         _norm = max(sum(abs(x')));
       } \
       else if (t == Inf) 
       {
         _norm = max(sum(abs(x)));
       } \
       else if (t == 2)
       {
         y  = rows 1 % 1;
         xt = x' * x;
         k  = 0;
         l  = 1/0;
         l1 = 0;
         while(abs(l-l1) >= eps & k < 100) 
         {
           z = xt * y;
           l1 = l;
           l = (y'*z)/(y'*y);
           m = 1/sqrt(sum(abs(z)^2));
           y = m * z;
           k = k + 1;
         }
         iter  = k;
         _norm = sqrt(l);
       } \
       else if (t == 2.1)
       {
         _norm = sqrt(max(eig(x'*x)[0:rows-1,0]));
       } \
       else if (t == "f")
       {
         _norm = sqrt(sum(diag(x'*x)));
       }
     } \
     else
     {
       if (t == Inf)
       {
         _norm = max(abs(x));
       } \
       else if (t == -Inf)
       {
         _norm = min(abs(x));
       } \
       else if (t <> 0)
       {
         _norm = sum(abs(x)^t)^(1/t);
       }
     }
   }
}

function jac(a, b)
import eps
export eps
{
  if (~exists("a"))
  {
    "Usage: jac(a[,b])";
  } \
  else
  {
    n = size(a)[0];
    if (~exists("b")) b = eye(size(a)[0]);
    if (~exists("eps")) eps = 1.0e-16;
    ev = jacob(a, b, eps);
    _jac = eigv;
    _jac[0:n-1,n] = ev';
  }
}

function triu(a)
{
  n = size(a)[0];
  m = size(a)[1];
  t = min(size(a));
  _triu = n m % 0;
  i = 0;
  for(i=0:t-1)
  {
    _triu[i,i:m-1] = a[i,i:m-1];
  }
}

function tril(a)
{
  n = size(a)[0];
  m = size(a)[1];
  t = min(size(a));
  _tril = n m % 0;
  for(i=0:t-1)
  {
    _tril[i:n-1,i] = a[i:n-1,i];
  }
}

function sort(x)
import Inf
{
  n = size(x)[0] * size(x)[1];
  _sort = x;
  if (n > 1) 
  {
    t = 1 n % x;
    v = 0:n-1;
    i = 0;
    while(i < n)
    {
      m = min(t);
      w = t[t==m];
      j = size(w)[1];
      y[i:i+j-1] = w;
      t[t==m] = Inf;
      i = i + j;
    }
    _sort = size(x) % y;
  }
}

function plot(y, x)
! PLOT curves from argument Y as function of argument X
! if argument X is not present use Y:s column index for 
! X axis.
!
! Curves plotted are rows of the argument matrix Y.
!
! Device can be selected by setting of variable gdevice, 
! x- and y axis lablebs by gxlabel and gylabel respectively
! and plot title by gtitle.
!
import gdevice,gxlabel,gylabel,gtitle
export gdevice
{
  ymax = max(max(y));
  ymin = min(min(y));

  np = size(y)[1];
  ny = np - 1;

  if (~exists("x"))
  {
    xmin = 1;
    xmax = np; 
    x=1:xmax;
  } else
  {
    xmin = min(x);
    xmax = max(x);
  }

  if ( ~exists("gdevice") ) askgdev
  gopen( gdevice );

  gwindow(0 1 0 1 0 1);   ! Set default window and viewport for axis drawing
  gviewport(0 1 0 1);

  gclear;                 ! Clear screen

  gcolor(1);              ! Use colorindex one

  gmove(0 0 0);           ! Plot frame around coordinate space
  gdraw(1 0 0);
  gdraw(1 1 0);
  gdraw(0 1 0);
  gdraw(0 0 0);

  gmove(0.1 0.9 0);       ! Axis vectors
  gdraw(0.1 0.1 0);
  gdraw(0.9 0.1 0);

  i = 0.9;                ! (X Axis) tick marks and values
  t = xmax;
  s = (xmax-xmin)/8;
  while(i > 0.1)
  {
    gmove(i 0.1 0);
    gdraw(i 0.075 0);
    str=sprintf("%.3g",t);
    gmove((i-size(str)[1]*0.013/2) 0.05 0);
    gtext(0.013 0, str);
    t=t-s;
    i=i-0.1;
  }

  i = 0.9;                ! (Y Axis) tick marks and values
  t = ymax;
  s = (ymax-ymin)/8;
  while(i > 0.1)
  {
    gmove(0.1 i 0);
    gdraw(0.075 i 0);
    str=sprintf("%.3g",t);
    gmove(0.05 (i-size(str)[1]*0.013/2) 0);
    gtext(0.013 90, str);
    t=t-s;
    i=i-0.1;
  }

  if (exists("gxlabel")) {
     gmove((0.5-size(gxlabel)[1]*0.019/2) 0.025 0);     ! X Axis label
     gtext(0.019 0, gxlabel);
  }

  if (exists("gylabel")) {
      gmove(0.025 (0.5-size(gylabel)[1]*0.019/2) 0);     ! Y Axis label
      gtext(0.019 90, gylabel);
  }

  if (~exists("gtitle")) gtitle="MATC PLOT";       ! Title
  gmove((0.5-size(gtitle)[1]*0.03/2) 0.925 0);
  gtext(0.03 0, gtitle);

  gwindow(xmin xmax ymin ymax 0 1); ! Set world coordinates
  gviewport(0.1 0.9 0.1 0.9);

  n = size(y)[0]; 
  p[0:ny,0]=x;
  p[0:ny,2]=0;
  for(i=0:n-1)                      ! Plot curves one at a time; each 
  {                                 ! with a new color
    p[0:ny,1]=y[i,0:ny];
    gcolor(i+1);
    gpolyline(np, p); 
  }

  gclose;
}

function contour( mat, scale, vscale )
!
! Plot equal value color contours from matrix mat. Axis numbering scales
! can be given with second parameter.
!
! Device can be selected by setting of variable gdevice, 
! x- and y axis lablebs by gxlabel and gylabel respectively
! and plot title by gtitle.
!
import gdevice,gxlabel,gylabel,gtitle
export gdevice
{
  if (~exists("scale"))
  {
    ymin = 1;
    ymax = size(mat)[0];
    xmin = 1;
    xmax = size(mat)[1];
  } else
  {
    xmin = scale[0];
    xmax = scale[1];
    ymin = scale[2];
    ymax = scale[3];
  }
 
  if ( ~exists("vscale") ) vscale=1

  if (~exists("gdevice")) askgdev
  gopen(gdevice);

  gwindow(0 1 0 1 0 1);   ! Set default window and viewport for axis drawing
  gviewport(0 1 0 vscale);

  gclear;                 ! Clear screen

  gcolor(1);              ! Use colorindex one

  gmove(0 0 0);           ! Plot frame around coordinate space
  gdraw(1 0 0);
  gdraw(1 1 0);
  gdraw(0 1 0);
  gdraw(0 0 0);

  gmove(0.1 0.9 0);       ! Axis vectors
  gdraw(0.1 0.1 0);
  gdraw(0.9 0.1 0);

  i = 0.9;                 ! (X Axis) tick marks and values
  t = xmax;
  s = (xmax-xmin)/8;
  while(i > 0.1)
  {
    gmove(i 0.1 0);
    gdraw(i 0.075 0);
    str=sprintf("%.3g",t);
    gmove((i-size(str)[1]*0.013/2) 0.05 0);
    gtext(0.013 0, str);
    t=t-s;
    i=i-0.1;
  }

  i = 0.9;                ! (Y Axis) tick marks and values
  t = ymax;
  s = (ymax-ymin)/8;
  while(i > 0.1)
  {
    gmove(0.1 i 0);
    gdraw(0.075 i 0);
    str=sprintf("%.3g",t);
    gmove(0.05 (i-size(str)[1]*0.013/2) 0);
    gtext(0.013 90, str);
    t=t-s;
    i=i-0.1;
  }

  if (exists("gxlabel")) {
      gmove((0.5-size(gxlabel)[1]*0.019/2) 0.025 0);     ! X Axis label
      gtext(0.019 0, gxlabel);
  }

  if (exists("gylabel")) {
      gmove(0.025 (0.5-size(gylabel)[1]*0.019/2) 0);     ! Y Axis label
      gtext(0.019 90, gylabel);
  }

  if (~exists("gtitle")) gtitle="MATC C3D/CONTOUR"; ! Title
  gmove((0.5-size(gtitle)[1]*0.03/2) 0.925 0);
  gtext(0.03 0, gtitle);

  gviewport(0.1 0.9 0.1 (vscale-0.1) );

  !
  ! array of pretty color selections
  !
  cols = 14 3 %   0   0 100 \
                  0  33 100 \
                  0  66 100 \
                  0 100 100 \
                  0  33   0 \
                  0  66   0 \
                  0 100   0 \
                100 100   0 \
                100 100   0 \
                100  66   0 \
                100  33   0 \
                100   0   0 \
                100   0 100 \
                100  33 100 \
                100  66 100;

  ! load colors
  for(i=2:15)
  {
    gdefcolor(i, cols[i-2,0 1 2]/100.0); 
  }

  gc3dlevels(13);       ! number of color levels
  gc3d(mat);            ! do it 

  gclose;
}

function mesh(mat, vp)
!
! Plot 3D mesh from matrix mat. Viewpoint can be 
! given with optional second parameter. Plot device is selected by
! imported variable gdevice, the title of the plot can be changed by
! variable gtitle.
!
import gdevice,gtitle
export gdevice
{
  if (~exists("gdevice")) askgdev
  gopen(gdevice);

  gwindow(0 1 0 1 0 1);   ! Set default window and viewport for axis drawing
  gviewport(0 1 0 1);

  gclear;                 ! Clear screen

  gcolor(1);              ! Use colorindex one

  gmove(0 0 0);           ! Plot frame around coordinate space
  gdraw(1 0 0);
  gdraw(1 1 0);
  gdraw(0 1 0);
  gdraw(0 0 0);

  if (~exists("gtitle")) gtitle="MATC C3D/MESH"; ! Title
  gmove((0.5-size(gtitle)[1]*0.03/2) 0.925 0);
  gtext(0.03 0, gtitle);

  gviewport(0.1 0.9 0.1 0.9);

  !  
  ! Set default or user given viewpoint 
  !
  if (exists("vp"))
  {
    gviewpoint(vp);
  } else
  {
    gviewpoint(1 3 % max(max(mat)));
  }

  gmove(0 0 0);
  gdraw(0 0 1);
  gmove(0 0 0);
  gdraw(0 1 0);
  gmove(0 0 0);
  gdraw(1 0 0);

  gc3dlevels(0);         ! number of color levels
  gc3d(mat);             ! do it

  gclose;
}

function cmesh(mat, vp)
!
! Plot 3D mesh and its equal value color contours from matrix mat.
! Viewpoint can be given with optional second parameter.
! Plot device is selected by imported variable gdevice, the 
! title of the plot can be changed by setting of gtitle.
!
import gdevice,gtitle
export gdevice
{
  if (~exists("gdevice")) askgdev;
  gopen(gdevice);

  gwindow(0 1 0 1 0 1);   ! Set default window and viewport for axis drawing
  gviewport(0 1 0 1);

  gclear;                 ! Clear screen

  gcolor(1);              ! Use colorindex one

  gmove(0 0 0);
  gdraw(1 0 0);
  gdraw(1 1 0);
  gdraw(0 1 0);
  gdraw(0 0 0);

  if (~exists("gtitle")) gtitle="MATC C3D/CMESH"; ! Title
  gmove((0.5-size(gtitle)[1]*0.03/2) 0.925 0);
  gtext(0.03 0, gtitle);

  gviewport(0.1 0.9 0.1 0.9);

  !
  ! array of pretty color selections
  !
  cols = 14 3 %   0   0 100 \
                  0  33 100 \
                  0  66 100 \
                  0 100 100 \
                  0  33   0 \
                  0  66   0 \
                  0 100   0 \
                100 100   0 \
                100 100   0 \
                100  66   0 \
                100  33   0 \
                100   0   0 \
                100   0 100 \
                100  66 100;

  ! load colors
  for(i=2:15)
  {
    gdefcolor(i, cols[i-2,0 1 2]/100.0); 
  }

  !  
  ! Set default or user given viewpoint 
  !
  if (exists("vp")) 
  {
    gviewpoint(vp); 
  } else
  {
    gviewpoint(1 3 % max(max(mat)));
  }

  gc3dlevels(13);      ! number of color levels
  gc3d(mat);           ! do it 

  gclose;
}

function printf(fmt,values)
{
  if (exists("values"))
  {
    fprintf(stdout,fmt,values);
  } else
  {
    fprintf(stdout,fmt);
  }
}

function scanf(fmt)
{
  _scanf = fscanf(stdin,fmt);
}

function gets(str)
{
  _gets = fgets(stdin);
}

function puts(str)
{
  fputs(stdout,str);
}

function tovector(a,b,c)
{
   n=0:size(a)[1]-1;
   _tovector(0,n)=a;
   _tovector(1,n)=b;
   if (exists("c"))
   {
     _tovector(2,n)=c;
   } else
   {
     _tovector(2,n)=0;
   }
}

NaN = ln(-1);
Inf = 1/0;
eps = 1.0e-15;
! pi  = 4 * atan(1);

format(6,1)
gdevice = 4

! source("./.mcinit");

fclose(stderr);