/****************************************************************/
/*								*/
/* xdmath.c:	stuff that depends on the			*/
/*		X Toolkit					*/
/*		and does data manipulation			*/
/*								*/
/*		Design: Walter Benzing 1994			*/
/*		NO RIGHTS RESERVED				*/
/*								*/
/****************************************************************/

/* standard includes */
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

/*includes for the X11 toolkit*/
#include <X11/Intrinsic.h>  /* Intrinsics Definitions */
#include <X11/StringDefs.h>
#include <X11/Shell.h>      /* Top Level Shell Widget */

/*my includes */
#include "xdim.h"

/*
 * Miror data set in X-direction
 */

void MA_MirrorX(Widget w, XtPointer client_data, XtPointer call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *view;
    int i, off;
    double *lp, *dp, *ep, val;
    unsigned char *clp, *cp, *cep, cval;

    XDim = (XDimInfo*)client_data;
    if(XDim->port >=0)
      {
       view = &((XDim->Views)[XDim->port]);
       view->changed |= DATA_AND_CHANGE;
       off = view->Width - 1;
       switch (view->dataType)
          {
           case MY_UCHAR:
	   if(view->charData == NULL)
             return;
	   clp = view->charData;
	   for(i = 0; i<view->Height; i++)
              {
               cp = clp;
               cep = cp+off;
               clp += view->Width;
               while(cep>cp)
        	  {
        	   cval = *cp;
        	   *cp++ = *cep;
        	   *cep = cval;
        	   cep--;
        	  }
              }
           break;
           default:
	   if(view->actData == NULL)
             return;
	   lp = view->actData;
	   for(i = 0; i<view->Height; i++)
              {
               dp = lp;
               ep = dp+off;
               lp += view->Width;
               while(ep>dp)
        	  {
        	   val = *dp;
        	   *dp++ = *ep;
        	   *ep = val;
        	   ep--;
        	  }
              }
	   break;
	  }
       XT_Redraw(((XDim->myWidgets).Views[XDim->port]).DArea,
	   (XtPointer) view, NULL);
      }
   }

/*
 * Miror data set in Y-direction
 */

void MA_MirrorY(Widget w, XtPointer client_data, XtPointer call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *view;
    int i, off, width;
    double *lp, *dp, *ep, val;
    unsigned char *clp, *cp, *cep, cval;

    XDim = (XDimInfo*)client_data;
    if(XDim->port >=0)
      {
       view = &((XDim->Views)[XDim->port]);
       view->changed |= DATA_AND_CHANGE;
       width = view->Width;
       off = (view->Height-1)*width;
       switch (view->dataType)
          {
           case MY_UCHAR:
	   if(view->charData == NULL)
             return;
	   clp = view->charData;
	   for(i = 0; i<width; i++)
              {
               cp = clp;
               clp ++;
               cep = cp+off;
               while(cep>cp)
        	  {
        	   cval = *cp;
        	   *cp = *cep;
        	   cp += width;
        	   *cep = cval;
        	   cep -= width;
        	  }
              }
           break;
           default:
	   if(view->actData == NULL)
             return;
	   lp = view->actData;
	   for(i = 0; i<width; i++)
              {
               dp = lp;
               lp ++;
               ep = dp+off;
               while(ep>dp)
        	  {
        	   val = *dp;
        	   *dp = *ep;
        	   dp += width;
        	   *ep = val;
        	   ep -= width;
        	  }
              }
           break;
          }
       if(w != NULL)
         XT_Redraw(((XDim->myWidgets).Views[XDim->port]).DArea,
	   (XtPointer) view, NULL);
      }
   }

/*
 * Negate data set
 */

void MA_Neg(Widget w, XtPointer client_data, XtPointer call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *view;
    int i,size;
    double *dp, val;

    XDim = (XDimInfo*)client_data;
    if(XDim->port >=0)
      {
       view = &((XDim->Views)[XDim->port]);
       if(view->dataType)
         WD_PopTypeCast(MY_DOUBLE, XDim->port, XDim);
       if(view->actData == NULL)
         return;
       view->changed |= DATA_AND_CHANGE;
       dp = view->actData;
       size = view->Width*view->Height;
       for(i = 0; i<size; i++)
          {
	   val = *dp;
	   *dp++ = -val;
          }
       val = view->Max;
       view->Max = -view->Min;
       view->Min = -val;
       XT_Redraw(((XDim->myWidgets).Views[XDim->port]).DArea,
	   (XtPointer) view, NULL);
      }
   }

/*
 * ABS
 */

void MA_Abs(Widget w, XtPointer client_data, XtPointer call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *view1, *view2;
    int i,size;
    double *dp1, *dp2, val;

    XDim = (XDimInfo*)client_data;
    if(XDim->port >=0)
      {
       view1 = &((XDim->Views)[XDim->port]);
       if(view1->dataType)
         WD_PopTypeCast(MY_DOUBLE, XDim->port, XDim);
       if(view1->actData == NULL)
         return;
       view2 = WD_CopyViewport("abs(", "", ")", view1, NULL, XDim);
       if(view2)
         {
	  dp1 = view1->actData;
	  dp2 = view2->actData;
	  size = view2->Width*view2->Height;
	  for(i = 0; i<size; i++)
             {
	      val = *dp1++;
	      *dp2++ = fabs(val);
             }
          CO_CalculateBases(view2);
	  XT_Redraw(((XDim->myWidgets).Views[XDim->port]).DArea,
	  (XtPointer) view2, NULL);
	 }
      }
   }

/*
 * Quadrate
 */

void MA_Quad(Widget w, XtPointer client_data, XtPointer call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *view1, *view2;
    int i,size;
    double *dp1, *dp2, val;

    XDim = (XDimInfo*)client_data;
    if(XDim->port >=0)
      {
       view1 = &((XDim->Views)[XDim->port]);
       if(view1->dataType)
         WD_PopTypeCast(MY_DOUBLE, XDim->port, XDim);
       if(view1->actData == NULL)
         return;
       view2 = WD_CopyViewport("(", "", ")^2", view1, NULL, XDim);
       if(view2)
         {
	  dp1 = view1->actData;
	  dp2 = view2->actData;
	  size = view2->Width*view2->Height;
	  for(i = 0; i<size; i++)
             {
	      val = *dp1++;
	      *dp2++ = val*val;
             }
          CO_CalculateBases(view2);
	  XT_Redraw(((XDim->myWidgets).Views[XDim->port]).DArea,
	  (XtPointer) view2, NULL);
	 }
      }
   }

/*
 * Square root
 */

void MA_Sqrt(Widget w, XtPointer client_data, XtPointer call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *view1, *view2;
    int i,size;
    double *dp1, *dp2, val;

    XDim = (XDimInfo*)client_data;
    if(XDim->port >=0)
      {
       view1 = &((XDim->Views)[XDim->port]);
       if(view1->dataType)
         WD_PopTypeCast(MY_DOUBLE, XDim->port, XDim);
       if(view1->actData == NULL)
         return;
       if(view1->Min < 0.0)
         {
          sprintf(XDim->error, "For sqrt the range must be x >= 0.0 !!");
          WD_XDimError(XDim, NULL);
          return;
         }
       view2 = WD_CopyViewport("sqrt(", "", ")", view1, NULL, XDim);
       if(view2)
         {
	  dp1 = view1->actData;
	  dp2 = view2->actData;
	  size = view2->Width*view2->Height;
	  for(i = 0; i<size; i++)
             {
	      val = *dp1++;
	      *dp2++ = sqrt(val);
             }
          CO_CalculateBases(view2);
	  XT_Redraw(((XDim->myWidgets).Views[XDim->port]).DArea,
	  (XtPointer) view2, NULL);
	 }
      }
   }

/*
 * log
 */

void MA_Log(Widget w, XtPointer client_data, XtPointer call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *view1, *view2;
    int i,size;
    double *dp1, *dp2, val;

    XDim = (XDimInfo*)client_data;
    if(XDim->port >=0)
      {
       view1 = &((XDim->Views)[XDim->port]);
       if(view1->dataType)
         WD_PopTypeCast(MY_DOUBLE, XDim->port, XDim);
       if(view1->actData == NULL)
         return;
       if(view1->Min <= 0.0)
         {
          sprintf(XDim->error, "For log the range must be x > 0.0 !!");
          WD_XDimError(XDim, NULL);
          return;
         }
       view2 = WD_CopyViewport("log(", "", ")", view1, NULL, XDim);
       if(view2)
         {
	  dp1 = view1->actData;
	  dp2 = view2->actData;
	  size = view2->Width*view2->Height;
	  for(i = 0; i<size; i++)
             {
	      val = *dp1++;
	      *dp2++ = log(val);
             }
          CO_CalculateBases(view2);
	  XT_Redraw(((XDim->myWidgets).Views[XDim->port]).DArea,
	  (XtPointer) view2, NULL);
	 }
      }
   }

/*
 * log10
 */

void MA_Log10(Widget w, XtPointer client_data, XtPointer call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *view1, *view2;
    int i,size;
    double *dp1, *dp2, val;

    XDim = (XDimInfo*)client_data;
    if(XDim->port >=0)
      {
       view1 = &((XDim->Views)[XDim->port]);
       if(view1->dataType)
         WD_PopTypeCast(MY_DOUBLE, XDim->port, XDim);
       if(view1->actData == NULL)
         return;
       if(view1->Min <= 0.0)
         {
          sprintf(XDim->error, "For log10 the range must be x > 0.0 !!");
          WD_XDimError(XDim, NULL);
          return;
         }
       view2 = WD_CopyViewport("log10(", "", ")", view1, NULL, XDim);
       if(view2)
         {
	  dp1 = view1->actData;
	  dp2 = view2->actData;
	  size = view2->Width*view2->Height;
	  for(i = 0; i<size; i++)
             {
	      val = *dp1++;
	      *dp2++ = log10(val);
             }
          CO_CalculateBases(view2);
	  XT_Redraw(((XDim->myWidgets).Views[XDim->port]).DArea,
	  (XtPointer) view2, NULL);
	 }
      }
   }

/*
 * exp
 */

void MA_Exp(Widget w, XtPointer client_data, XtPointer call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *view1, *view2;
    int i,size;
    double *dp1, *dp2, val;

    XDim = (XDimInfo*)client_data;
    if(XDim->port >=0)
      {
       view1 = &((XDim->Views)[XDim->port]);
       if(view1->dataType)
         WD_PopTypeCast(MY_DOUBLE, XDim->port, XDim);
       if(view1->actData == NULL)
         return;
       view2 = WD_CopyViewport("exp(", "", ")", view1, NULL, XDim);
       if(view2)
         {
	  dp1 = view1->actData;
	  dp2 = view2->actData;
	  size = view2->Width*view2->Height;
	  for(i = 0; i<size; i++)
             {
	      val = *dp1++;
	      *dp2++ = exp(val);
             }
          CO_CalculateBases(view2);
	  XT_Redraw(((XDim->myWidgets).Views[XDim->port]).DArea,
	  (XtPointer) view2, NULL);
	 }
      }
   }

/*
 * 1/x
 */

void MA_1x(Widget w, XtPointer client_data, XtPointer call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *view1, *view2;
    int i,size;
    double *dp1, *dp2, val;

    XDim = (XDimInfo*)client_data;
    if(XDim->port >=0)
      {
       view1 = &((XDim->Views)[XDim->port]);
       if(view1->dataType)
         WD_PopTypeCast(MY_DOUBLE, XDim->port, XDim);
       if(view1->actData == NULL)
         return;
       size = view1->Width*view1->Height;
       dp1 = view1->actData;
       for(i = 0; i<size; i++)
          {
           if(*dp1++ == 0.0)
             break;
          }
       if(i<size)
         {
          strcpy(XDim->error, "For 1/x all x must be != 0.0 !!");
          WD_XDimError(XDim, NULL);
          return;
         }
       view2 = WD_CopyViewport("1/(", "", ")", view1, NULL, XDim);
       if(view2)
         {
	  dp1 = view1->actData;
	  dp2 = view2->actData;
	  for(i = 0; i<size; i++)
             {
	      val = *dp1++;
	      *dp2++ = 1/val;
             }
          CO_CalculateBases(view2);
	  XT_Redraw(((XDim->myWidgets).Views[XDim->port]).DArea,
	  (XtPointer) view2, NULL);
	 }
      }
   }

/*
 * Scale data set
 */

void MA_Scale(int num, void *arg, void *call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *view1, *view2;
    char *value, start[256];
    int i,size;
    double fact, *dp1, *dp2, val;

    XDim = (XDimInfo*)call_data;
    value = (char *) arg;
    if(!sscanf(value, "%le", &fact))
      {
       strcpy(XDim->error, "Scale factor must be given !!");
       WD_XDimError(XDim, NULL);
       return;
      }
    if(XDim->port >=0)
      {
       view1 = &((XDim->Views)[XDim->port]);
       if(view1->dataType)
         WD_PopTypeCast(MY_DOUBLE, XDim->port, XDim);
       if(view1->actData == NULL)
         return;
       sprintf(start, "%s*(", value);
       view2 = WD_CopyViewport(start, "", ")", view1, NULL, XDim);
       if(view2)
         {
	  dp1 = view1->actData;
	  dp2 = view2->actData;
	  size = view2->Width*view2->Height;
	  for(i = 0; i<size; i++)
             {
	      val = *dp1++;
	      *dp2++ = fact*val;
             }
          CO_CalculateBases(view2);
	  XT_Redraw(((XDim->myWidgets).Views[XDim->port]).DArea,
	  (XtPointer) view2, NULL);
	 }
      }
   }

/*
 * get scale value
 */

void MA_GetScale(Widget w, XtPointer client_data, XtPointer call_data)
   {
    XDimInfo *XDim;

    XDim = (XDimInfo*)client_data;
    WD_PopInput("Scale Factor :", "1.0", 5, 10, MA_Scale, XDim);
   }

/*
 * Add an offset
 */

void MA_Offset(int num, void *arg, void *call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *view1, *view2;
    char *value;
    myString start;
    int i,size;
    double fact, *dp1, *dp2, val;

    XDim = (XDimInfo*)call_data;
    value = (char *) arg;
    if(!sscanf(value, "%le", &fact))
      {
       strcpy(XDim->error, "Offset must be given !!");
       WD_XDimError(XDim, NULL);
       return;
      }
    if(XDim->port >=0 && fact != 0.0)
      {
       view1 = &((XDim->Views)[XDim->port]);
       if(view1->dataType)
         WD_PopTypeCast(MY_DOUBLE, XDim->port, XDim);
       if(view1->actData == NULL)
         return;
       if(fact > 0)
         sprintf(start, ")+%s", value);
       else
         sprintf(start, ")%s", value);
       view2 = WD_CopyViewport("(", "", start, view1, NULL, XDim);
       if(view2)
         {
	  dp1 = view1->actData;
	  dp2 = view2->actData;
	  size = view2->Width*view2->Height;
	  for(i = 0; i<size; i++)
             {
	      val = *dp1++;
	      *dp2++ = fact+val;
             }
          CO_CalculateBases(view2);
	  XT_Redraw(((XDim->myWidgets).Views[XDim->port]).DArea,
	  (XtPointer) view2, NULL);
	 }
      }
   }

/*
 * get Offset value
 */

void MA_GetOffset(Widget w, XtPointer client_data, XtPointer call_data)
   {
    XDimInfo *XDim;

    XDim = (XDimInfo*)client_data;
    WD_PopInput("Offset value:", "0.0", 5, 10, MA_Offset, XDim);
   }


/*
 * edit value
 */

void MA_okEditVal(int num, void *arg, void *call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *view;
    char *value;
    int size, diff;
    double fact;

    XDim = (XDimInfo*)call_data;
    value = (char *) arg;
    if(!sscanf(value, "%le", &fact))
      {
       strcpy(XDim->error, "No value given !!");
       WD_XDimError(XDim, NULL);
       return;
      }
    if(XDim->port >=0 && XDim->aView.edPort >= 0)
      {
       view = &((XDim->Views)[XDim->aView.edPort]);
       diff = XDim->aView.edit - view->actData;
       size = view->Width * view->Height;
       if(view->actData == NULL || view->pict == NULL
		|| diff < 0 || diff >= size)
         return;
       *(XDim->aView.edit) = fact;
       view->changed |= DATA_AND_CHANGE;
       if(fact < view->Min || fact > view->Max)
         CO_CalculateBases(view);
       XT_Redraw(((XDim->myWidgets).Views[view->port]).DArea,
       (XtPointer) view, NULL);
      }
   }

/*
 * interpolate data field
 */

void MA_Interpolate(int newWidth, int newHeight, ViewportInfo *view1,
	XDimInfo *XDim)
   {
    ViewportInfo*view2;
    int i, j, port;
    double *dp, x, y, factX, factY;

    if(view1->dataType)
      WD_PopTypeCast(MY_DOUBLE, XDim->port, XDim);
    if(XDim->port >= 0 && view1->actData &&
       ((newWidth != view1->Width) || (newHeight != view1->Height)))
      {
       port = XDim->port;
       WD_CreateViewport((XDim->myWidgets).NewView, XDim, NULL);
       if(port != XDim->port)
	 {
	  if(CO_NewData(newWidth, newHeight, MY_DOUBLE, XDim))
	    {
	     WD_XDimError(XDim, NULL);
	     XT_DestroyViewport(NULL,&(XDim->port), NULL); 
	    }
	  else
	    {
             view2 = &((XDim->Views)[XDim->port]);
	     strcpy(view2->DimX, view1->DimX);
	     strcpy(view2->DimY, view1->DimY);
	     strcpy(view2->DimZ, view1->DimZ);
	     if(newWidth > 1)
	       factX = ((double)(view1->Width-1))/(newWidth-1);
	     else
	       factX = 1.0;
	     if(newHeight > 1)
	       factY = ((double)(view1->Height-1))/(newHeight-1);
	     else
	       factY = 1.0;
	     view2->FactX = factX*view1->FactX;
	     view2->FactY = factY*view1->FactY;
	     view2->OffX = view1->OffX;
	     view2->OffY = view1->OffY;
	     dp = view2->actData;
	     y = 0.0;
	     for(i=0; i<newHeight; i++)
	        {
	         x = 0.0;
	         for(j=0; j<newWidth; j++)
	            {
		     *dp++ = CO_Inter(x,y,view1);
		     x += factX;
		    }
	         y += factY;
		}
             CO_CalculateBases(view2);
	     XT_Redraw(((XDim->myWidgets).Views[XDim->port]).DArea,
	     (XtPointer) view2, NULL);
	    }
	 }       
      }
   }

/*
 * interpolate data field with spline
 */

void MA_Spline(int newWidth, int newHeight, ViewportInfo *view1,
	XDimInfo *XDim)
   {
    ViewportInfo*view2;
    int i, j, port, size;
    double *dp, *dp1, *dp2, *sx, *sy, *zp, x, y, factX, factY;

    if(view1->dataType)
      WD_PopTypeCast(MY_DOUBLE, XDim->port, XDim);
    if(XDim->port >= 0 && view1->actData &&
       ((newWidth != view1->Width) || (newHeight != view1->Height)))
      {
       port = XDim->port;
       WD_CreateViewport((XDim->myWidgets).NewView, XDim, NULL);
       if(port != XDim->port)
	 {
	  if(CO_NewData(newWidth, newHeight, MY_DOUBLE, XDim))
	    {
	     WD_XDimError(XDim, NULL);
	     XT_DestroyViewport(NULL,&(XDim->port), NULL); 
	    }
	  else
	    {
             view2 = &((XDim->Views)[XDim->port]);
	     strcpy(view2->DimX, view1->DimX);
	     strcpy(view2->DimY, view1->DimY);
	     strcpy(view2->DimZ, view1->DimZ);
	     if(newWidth > 1)
	       factX = ((double)(view1->Width-1))/(newWidth-1);
	     else
	       factX = 1.0;
	     if(newHeight > 1)
	       factY = ((double)(view1->Height-1))/(newHeight-1);
	     else
	       factY = 1.0;
	     view2->FactX = factX*view1->FactX;
	     view2->FactY = factY*view1->FactY;
	     view2->OffX = view1->OffX;
	     view2->OffY = view1->OffY;
	     size = view1->Width*view1->Height*sizeof(double);
	     sx = (double*)malloc(size);
	     size = newWidth*view1->Height*sizeof(double);
	     sy = (double*)malloc(size);
	     zp = (double*)malloc(size);
	     if(sx==NULL || sy == NULL || zp == NULL)
	       {
	        sprintf(XDim->error, "Couldn't allocate %d number of bytes",
			2*size+view1->Width*view1->Height);
		WD_XDimError(XDim, NULL);
	        XT_DestroyViewport(NULL,&(XDim->port), NULL); 
	       }
	     else
	       {
		dp = view1->actData;
		dp1 = sx;
		for(i=0; i<view1->Height; i++)
	           {
	            if(CO_StartSpline(view1->Width, 1.0, dp1, dp, XDim))
	               {
	                WD_XDimError(XDim, NULL);
	                XT_DestroyViewport(NULL, &(XDim->port), NULL);
			free(sx);
			free(sy);
			free(zp);
			return;
		       }
	            dp1 += view1->Width;
	            dp += view1->Width;
		   }
		dp = zp;
	        x = 0.0;
		for(i=0; i<newWidth; i++)
	           {
	            dp1 = sx;
	            dp2 = view1->actData;
	            for(j=0; j<view1->Height; j++)
	               {
	                *dp++ = CO_Spline(x, 1.0, view1->Width, dp1, dp2);
	        	dp1 += view1->Width;
	        	dp2 += view1->Width;
	               }
	            x += factX;
		   }
		dp = zp;
		dp1 = sy;
		for(i=0; i<newWidth; i++)
	           {
	            if(CO_StartSpline(view1->Height, 1.0, dp1, dp, XDim))
	               {
	                WD_XDimError(XDim, NULL);
	                XT_DestroyViewport(NULL, &(XDim->port), NULL);
			free(sx);
			free(sy);
			free(zp);
			return;
		       }
	            dp1 += view1->Height;
	            dp += view1->Height;
		   }
		dp = view2->actData;
	        y = 0.0;
		for(i=0; i<newHeight; i++)
	           {
	            dp1 = sy;
	            dp2 = zp;
	            for(j=0; j<newWidth; j++)
	               {
	                *dp++ = CO_Spline(y, 1.0, view1->Height, dp1, dp2);
	        	dp1 += view1->Height;
	        	dp2 += view1->Height;
	               }
	            y += factY;
		   }
		free(sx);
		free(sy);
		free(zp);
        	CO_CalculateBases(view2);
		XT_Redraw(((XDim->myWidgets).Views[XDim->port]).DArea,
		(XtPointer) view2, NULL);
	       }
	    }
	 }       
      }
   }


/*
 * Now operations with fields
 */

/*
 * create new data field
 */

void MA_New(int num, void *arg, void *call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *view1, *view2;
    myString *vals;
    int width, height, size, port, i;
    double *dp;

    XDim = (XDimInfo*)call_data;
    vals = (myString *)arg;
    width =0;
    height = 0;
    sscanf(vals[0], "%d", &width);
    sscanf(vals[1], "%d", &height);
    if(width <= 0 || height <= 0)
      {
       strcpy(XDim->error, "Width and Height of the new field must be > 0");
       WD_XDimError(XDim, NULL);
      }
    else
      {
       port = XDim->port;
       WD_CreateViewport((XDim->myWidgets).NewView, XDim, NULL);
       if(port != XDim->port)
	 {
	  if(CO_NewData(width, height, MY_DOUBLE, XDim))
	    {
	     WD_XDimError(XDim, NULL);
	     XT_DestroyViewport(NULL,&(XDim->port), NULL); 
	    }
	  else
	    {
             view1 = &((XDim->Views)[port]);
             view2 = &((XDim->Views)[XDim->port]);
	     strcpy(view2->DimX, view1->DimX);
	     strcpy(view2->DimY, view1->DimY);
	     strcpy(view2->DimZ, view1->DimZ);
	     view2->FactX = view1->FactX;
	     view2->FactY = view1->FactY;
	     view2->OffX = view1->OffX;
	     view2->OffY = view1->OffY;
	     dp = view2->actData;
	     size = width * height;
	     for(i=0; i<size; i++)
	        *dp++ = 0.0;
	     CO_CalculateBases(view2);
	     XT_Redraw(((XDim->myWidgets).Views[XDim->port]).DArea,
	     (XtPointer) view2, NULL);
	    }
	 }       
      }
   }

/*
 * start new field
 */

void MA_StartNew(Widget w, XtPointer client_data, XtPointer call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *view;
    char width[10], height[10];

    
    XDim = (XDimInfo*)client_data;
    if(XDim->port >= 0)
      {
       view = &((XDim->Views)[XDim->port]);
       sprintf(width, "%d", view->Width);
       sprintf(height, "%d", view->Height);
       WD_PopInput2("Field Width:", width, 4, 4, "Field Height:", height, 4, 4,
	MA_New, XDim);
      }
   }

/*
 * add two data fields
 */

void MA_Add(int num, void *arg, void *call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *view1, *view2, *view3;
    int i, j, *ip, x0, y0, x1, y1;
    double *lp2, *dp1, *dp2, *dp3;

    XDim = (XDimInfo*)call_data;
    view1 = &((XDim->Views)[XDim->oper1]);
    view2 = &((XDim->Views)[XDim->oper2]);
    if(num)
      {
       ip = arg;
       x0 = ip[0];
       y0 = ip[1];
      }
    else
      {
       x0 = 0;
       y0 = 0;
      }
    view3 = WD_CopyViewport("(", ")+(", ")", view1, view2, XDim);
    if(view3)
      {
       dp1 = view1->actData;
       dp3 = view3->actData;
       y0 = -y0;
       x0 = -x0;
       lp2 = view2->actData+y0*view2->Width+x0;
       y1 = y0;
       for(i=0; i<view1->Height; i++)
          {
           dp2 = lp2;
           lp2 += view2->Width;
           x1 = x0;
           for(j=0; j < view1->Width; j++)
              {
               if(y1>=0 && y1 < view2->Height && x1 >=0 && x1<view2->Width)
                 *dp3++ = *dp1++ + *dp2;
               else
                 *dp3++ = *dp1++;
               dp2++;
	       x1++;               
              }
           y1++;
          }
       CO_CalculateBases(view3);
       XT_Redraw(((XDim->myWidgets).Views[XDim->port]).DArea,
       (XtPointer) view3, NULL);
      }
   }

/*
 * start add
 */
void MA_StartAdd(Widget w, XtPointer client_data, XtPointer call_data)
   {
    XDimInfo *XDim;

    XDim = (XDimInfo*)client_data;
    WD_PopMath("Add to Field", "Field", XDim->port, MA_Add, XDim);
   }

/*
 * Subtract two data fields
 */

void MA_Sub(int num, void *arg, void *call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *view1, *view2, *view3;
    int i, j, *ip, x0, y0, x1, y1;
    double *lp2, *dp1, *dp2, *dp3;

    XDim = (XDimInfo*)call_data;
    view1 = &((XDim->Views)[XDim->oper1]);
    view2 = &((XDim->Views)[XDim->oper2]);
    if(num)
      {
       ip = arg;
       x0 = ip[0];
       y0 = ip[1];
      }
    else
      {
       x0 = 0;
       y0 = 0;
      }
    view3 = WD_CopyViewport("(", ")-(", ")", view1, view2, XDim);
    if(view3)
      {
       dp1 = view1->actData;
       dp3 = view3->actData;
       y0 = -y0;
       x0 = -x0;
       lp2 = view2->actData+y0*view2->Width+x0;
       y1 = y0;
       for(i=0; i<view1->Height; i++)
          {
           dp2 = lp2;
           lp2 += view2->Width;
           x1 = x0;
           for(j=0; j < view1->Width; j++)
              {
               if(y1>=0 && y1 < view2->Height && x1 >=0 && x1<view2->Width)
                 *dp3++ = *dp1++ - *dp2;
               else
                 *dp3++ = *dp1++;
               dp2++;
	       x1++;               
              }
           y1++;
          }
       CO_CalculateBases(view3);
       XT_Redraw(((XDim->myWidgets).Views[XDim->port]).DArea,
       (XtPointer) view3, NULL);
      }
   }

/*
 * start add
 */

void MA_StartSub(Widget w, XtPointer client_data, XtPointer call_data)
   {
    XDimInfo *XDim;

    XDim = (XDimInfo*)client_data;
    WD_PopMath("Subtract from Field", "Field", XDim->port, MA_Sub, XDim);
   }

/*
 * Multiply two data fields
 */

void MA_Mul(int num, void *arg, void *call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *view1, *view2, *view3;
    int i, j, *ip, x0, y0, x1, y1;
    double *lp2, *dp1, *dp2, *dp3;

    XDim = (XDimInfo*)call_data;
    view1 = &((XDim->Views)[XDim->oper1]);
    view2 = &((XDim->Views)[XDim->oper2]);
    if(num)
      {
       ip = arg;
       x0 = ip[0];
       y0 = ip[1];
      }
    else
      {
       x0 = 0;
       y0 = 0;
      }
    view3 = WD_CopyViewport("(", ")*(", ")", view1, view2, XDim);
    if(view3)
      {
       dp1 = view1->actData;
       dp3 = view3->actData;
       y0 = -y0;
       x0 = -x0;
       lp2 = view2->actData+y0*view2->Width+x0;
       y1 = y0;
       for(i=0; i<view1->Height; i++)
          {
           dp2 = lp2;
           lp2 += view2->Width;
           x1 = x0;
           for(j=0; j < view1->Width; j++)
              {
               if(y1>=0 && y1 < view2->Height && x1 >=0 && x1<view2->Width)
                 *dp3++ = *dp1++ * *dp2;
               else
                 *dp3++ = *dp1++;
               dp2++;
	       x1++;               
              }
           y1++;
          }
       CO_CalculateBases(view3);
       XT_Redraw(((XDim->myWidgets).Views[XDim->port]).DArea,
       (XtPointer) view3, NULL);
      }
   }

/*
 * start Multiplication
 */

void MA_StartMul(Widget w, XtPointer client_data, XtPointer call_data)
   {
    XDimInfo *XDim;

    XDim = (XDimInfo*)client_data;
    WD_PopMath("Multiply Field", "by", XDim->port, MA_Mul, XDim);
   }

/*
 * Divide two data fields
 */

void MA_Div(int num, void *arg, void *call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *view1, *view2, *view3;
    int i, j, *ip, x0, y0, x1, y1;
    double *lp2, *dp1, *dp2, *dp3;
    Boolean Error = False;

    XDim = (XDimInfo*)call_data;
    view1 = &((XDim->Views)[XDim->oper1]);
    view2 = &((XDim->Views)[XDim->oper2]);
    if(num)
      {
       ip = arg;
       x0 = -ip[0];
       y0 = -ip[1];
      }
    else
      {
       x0 = 0;
       y0 = 0;
      }
    lp2 = view2->actData+y0*view2->Width+x0;
    y1 = y0;
    for(i=0; i<view1->Height; i++)
       {
        dp2 = lp2;
        lp2 += view2->Width;
        x1 = x0;
        for(j=0; j < view1->Width; j++)
           {
            if(y1>=0 && y1 < view2->Height && x1 >=0 && x1<view2->Width)
              if(*dp2++ == 0.0)
                {
                 Error = True;
                 i = view1->Height;
                 break;
                }
            else
	      dp2++;
	    x1++;               
           }
        y1++;
       }
    if(Error)
      {
       strcpy(XDim->error,
	     "To divide all selected field elements must be != 0");
       WD_XDimError(XDim, NULL);
       view3 = NULL;
      }
    else
      view3 = WD_CopyViewport("(", ")/(", ")", view1, view2, XDim);
    if(view3)
      {
       dp1 = view1->actData;
       dp3 = view3->actData;
       lp2 = view2->actData+y0*view2->Width+x0;
       y1 = y0;
       for(i=0; i<view1->Height; i++)
          {
           dp2 = lp2;
           lp2 += view2->Width;
           x1 = x0;
           for(j=0; j < view1->Width; j++)
              {
               if(y1>=0 && y1 < view2->Height && x1 >=0 && x1<view2->Width)
                 *dp3++ = *dp1++ / *dp2;
               else
                 *dp3++ = *dp1++;
               dp2++;
	       x1++;               
              }
           y1++;
          }
       CO_CalculateBases(view3);
       XT_Redraw(((XDim->myWidgets).Views[XDim->port]).DArea,
       (XtPointer) view3, NULL);
      }
   }

/*
 * start division
 */

void MA_StartDiv(Widget w, XtPointer client_data, XtPointer call_data)
   {
    XDimInfo *XDim;

    XDim = (XDimInfo*)client_data;
    WD_PopMath("Divide Field", "by", XDim->port, MA_Div, XDim);
   }
/*
 * Complex Math
 */

/*
 * XY->ABS/PHI
 */

void MA_AbsPhi(Widget w, XtPointer client_data, XtPointer call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *vreal, *vimag, *vabs, *vphi;
    int real, imag, abs, phi, i,size;
    double *dr, *di, *da, *dp, re, im;

    XDim = (XDimInfo*)client_data;
    if(XDim->port >=0)
      {
       vreal = &((XDim->Views)[XDim->port]);
       real = FI_GetLink("real", vreal, XDim);
       imag = FI_GetLink("imag", vreal, XDim);
       if(real < 0 || imag < 0)
         {
          strcpy(XDim->mess, "Link to REAL or IMAGinary part not found");
          WD_ViewportMessage(vreal, XDim);
	  return;
         }
       vreal = &((XDim->Views)[real]);
       if(vreal->dataType)
         {
          WD_PopTypeCast(MY_DOUBLE, real, XDim);
          return;
         }
       vimag = &((XDim->Views)[imag]);
       if(vimag->dataType)
         {
          WD_PopTypeCast(MY_DOUBLE, imag, XDim);
          return;
         }
       if(vreal->actData == NULL || vimag->actData == NULL)
         {
          strcpy(XDim->mess, "Couldn't process zero data field !");
          WD_ViewportMessage(vreal, XDim);
          WD_ViewportMessage(vimag, XDim);
          return;
         }
       if(vreal->Width != vimag->Width || vreal->Height != vimag->Height)
         {
          strcpy(XDim->mess, "Size of real and imag part must be the same !");
          WD_ViewportMessage(vreal, XDim);
          WD_ViewportMessage(vimag, XDim);
          return;
         }
       vabs = WD_CopyViewport("Abs(", "+i", ")", vreal, vimag, XDim);
       vphi = WD_CopyViewport("Phi(", "+i", ")", vreal, vimag, XDim);
       if(vabs != NULL && vphi != NULL)
         {
          abs = vabs->port;
          phi = vphi->port;
	  strcpy(vphi->DimZ, "rad");
	  dr = vreal->actData;
	  di = vimag->actData;
	  da = vabs->actData;
	  dp = vphi->actData;
	  size = vreal->Width*vreal->Height;
	  for(i = 0; i<size; i++)
             {
	      re = *dr++;
	      im = *di++;
	      *da++ = sqrt(re*re+im*im);
	      *dp++ = atan2(im,re);
             }
          CO_CalculateBases(vabs);
          CO_CalculateBases(vphi);
          FI_AddLink(abs, phi, "abs", "phi", XDim);
	  XT_Redraw(((XDim->myWidgets).Views[abs]).DArea,
	  (XtPointer) vabs, NULL);
	  XT_Redraw(((XDim->myWidgets).Views[phi]).DArea,
	  (XtPointer) vphi, NULL);
	 }
      }
   }

/*
 * ABS/PHI -> XY
 */

void MA_cxy(Widget w, XtPointer client_data, XtPointer call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *vreal, *vimag, *vabs, *vphi;
    int real, imag, abs, phi, i,size;
    double *dr, *di, *da, *dp, absVal, phiVal;

    XDim = (XDimInfo*)client_data;
    if(XDim->port >=0)
      {
       vreal = &((XDim->Views)[XDim->port]);
       abs = FI_GetLink("abs", vreal, XDim);
       phi = FI_GetLink("phi", vreal, XDim);
       if(abs < 0 || phi < 0)
         {
          strcpy(XDim->mess, "Link to ABS or PHI not found");
          WD_ViewportMessage(vreal, XDim);
	  return;
         }
       vabs = &((XDim->Views)[abs]);
       if(vabs->dataType)
         {
          WD_PopTypeCast(MY_DOUBLE, abs, XDim);
          return;
         }
       vphi = &((XDim->Views)[phi]);
       if(vphi->dataType)
         {
          WD_PopTypeCast(MY_DOUBLE, phi, XDim);
          return;
         }
       if(vabs->actData == NULL || vphi->actData == NULL)
         {
          strcpy(XDim->mess, "Couldn't process zero data field !");
          WD_ViewportMessage(vabs, XDim);
          WD_ViewportMessage(vphi, XDim);
          return;
         }
       if(vabs->Width != vphi->Width || vabs->Height != vphi->Height)
         {
          strcpy(XDim->mess, "Size of real and imag part must be the same !");
          WD_ViewportMessage(vabs, XDim);
          WD_ViewportMessage(vphi, XDim);
          return;
         }
       vreal = WD_CopyViewport("real(", "*e^i", ")", vabs, vphi, XDim);
       vimag = WD_CopyViewport("imag(", "*e^i", ")", vabs, vphi, XDim);
       if(vreal != NULL && vimag != NULL)
         {
          real = vreal->port;
          imag = vimag->port;
	  dr = vreal->actData;
	  di = vimag->actData;
	  da = vabs->actData;
	  dp = vphi->actData;
	  size = vabs->Width*vabs->Height;
	  for(i = 0; i<size; i++)
             {
	      absVal = *da++;
	      phiVal = *dp++;
	      if(phiVal >0)
	        phiVal *= 1.0;
	      *di++ = absVal*sin(phiVal);
	      *dr++ = absVal*cos(phiVal);
             }
          CO_CalculateBases(vreal);
          CO_CalculateBases(vimag);
          FI_AddLink(real, imag, "real", "imag", XDim);
	  XT_Redraw(((XDim->myWidgets).Views[real]).DArea,
	  (XtPointer) vreal, NULL);
	  XT_Redraw(((XDim->myWidgets).Views[imag]).DArea,
	  (XtPointer) vimag, NULL);
	 }
      }
   }

/*
 * add complex fields
 */

void MA_CAdd(int num, void *arg, void *call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *vreal1, *vimag1, *vreal2, *vimag2, *vreal3, *vimag3;
    int i, j, *ip, x0, y0, x1, y1, off;
    int real, imag;
    double *dr1, *dr2, *dr3, *di1, *di2, *di3;
    double *li2, *lr2;

    XDim = (XDimInfo*)call_data;
    vreal2 = &((XDim->Views)[XDim->oper2]);
    real = FI_GetLink("real", vreal2, XDim);
    imag = FI_GetLink("imag", vreal2, XDim);
    if(real < 0 || imag < 0)
      {
       strcpy(XDim->mess, "Link to REAL or IMAGinary part not found");
       WD_ViewportMessage(vreal2, XDim);
       return;
      }
    vreal2 = &((XDim->Views)[real]);
    if(vreal2->dataType)
      {
       WD_PopTypeCast(MY_DOUBLE, real, XDim);
       return;
      }
    vimag2 = &((XDim->Views)[imag]);
    if(vimag2->dataType)
      {
       WD_PopTypeCast(MY_DOUBLE, imag, XDim);
       return;
      }
    vreal1 = &((XDim->Views)[XDim->oper1]);
    real = FI_GetLink("real", vreal1, XDim);
    imag = FI_GetLink("imag", vreal1, XDim);
    vreal1 = &((XDim->Views)[real]);
    vimag1 = &((XDim->Views)[imag]);
    if(vreal1->actData == NULL || vimag1->actData == NULL)
      {
       strcpy(XDim->mess, "Couldn't process zero data field !");
       WD_ViewportMessage(vreal1, XDim);
       WD_ViewportMessage(vimag1, XDim);
       return;
      }
    if(vreal1->Width != vimag1->Width || vreal1->Height != vimag1->Height)
      {
       strcpy(XDim->mess, "Size of real and imag part must be the same !");
       WD_ViewportMessage(vreal1, XDim);
       WD_ViewportMessage(vimag1, XDim);
       return;
      }
    if(num)
      {
       ip = arg;
       x0 = -ip[0];
       y0 = -ip[1];
      }
    else
      {
       x0 = 0;
       y0 = 0;
      }
    vreal3 = WD_CopyViewport("(", "+", ")", vreal1, vreal2, XDim);
    vimag3 = WD_CopyViewport("(", "+", ")", vimag1, vimag2, XDim);
    if(vreal3 != NULL && vimag3 != NULL)
      {
       dr1 = vreal1->actData;
       di1 = vimag1->actData;
       dr3 = vreal3->actData;
       di3 = vimag3->actData;
       off = y0*vreal2->Width+x0;
       lr2 = vreal2->actData + off;
       li2 = vimag2->actData + off;
       y1 = y0;
       for(i=0; i<vreal1->Height; i++)
          {
           dr2 = lr2;
           di2 = li2;
           lr2 += vreal2->Width;
           li2 += vimag2->Width;
           x1 = x0;
           for(j=0; j < vreal1->Width; j++)
              {
               if(y1>=0 && y1 < vreal2->Height && x1 >=0 && x1<vreal2->Width)
                 {
                  *dr3++ = *dr1++ + *dr2++;
                  *di3++ = *di1++ + *di2++;
                 }
               else
                 {
                  *dr3++ = *dr1++;
                  *di3++ = *di1++;
        	  dr2++;
        	  di2++;
                 }
	       x1++;               
              }
           y1++;
          }
       CO_CalculateBases(vreal3);
       CO_CalculateBases(vimag3);
       FI_AddLink(vreal3->port, vimag3->port, "real", "imag", XDim);
       XT_Redraw(((XDim->myWidgets).Views[vreal3->port]).DArea,
       (XtPointer) vreal3, NULL);
       XT_Redraw(((XDim->myWidgets).Views[vimag3->port]).DArea,
       (XtPointer) vimag3, NULL);
      }
   }

/*
 * start complex add
 */

void MA_StartCAdd(Widget w, XtPointer client_data, XtPointer call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *vreal, *vimag;
    int real, imag;

    XDim = (XDimInfo*)client_data;
    vreal = &((XDim->Views)[XDim->port]);
    real = FI_GetLink("real", vreal, XDim);
    imag = FI_GetLink("imag", vreal, XDim);
    if(real < 0 || imag < 0)
      {
       strcpy(XDim->mess, "Link to REAL or IMAGinary part not found");
       WD_ViewportMessage(vreal, XDim);
       return;
      }
    vreal = &((XDim->Views)[real]);
    if(vreal->dataType)
      {
       WD_PopTypeCast(MY_DOUBLE, real, XDim);
       return;
      }
    vimag = &((XDim->Views)[imag]);
    if(vimag->dataType)
      {
       WD_PopTypeCast(MY_DOUBLE, imag, XDim);
       return;
      }
    if(vreal->actData == NULL || vimag->actData == NULL)
      {
       strcpy(XDim->mess, "Couldn't process zero data field !");
       WD_ViewportMessage(vreal, XDim);
       WD_ViewportMessage(vimag, XDim);
       return;
      }
    if(vreal->Width != vimag->Width || vreal->Height != vimag->Height)
      {
       strcpy(XDim->mess, "Size of real and imag part must be the same !");
       WD_ViewportMessage(vreal, XDim);
       WD_ViewportMessage(vimag, XDim);
       return;
      }
    WD_PopMath("Add to Complex Field", "Complex Field", real, MA_CAdd, XDim);
   }

/*
 * sub complex fields
 */

void MA_CSub(int num, void *arg, void *call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *vreal1, *vimag1, *vreal2, *vimag2, *vreal3, *vimag3;
    int i, j, *ip, x0, y0, x1, y1, off;
    int real, imag;
    double *dr1, *dr2, *dr3, *di1, *di2, *di3;
    double *li2, *lr2;

    XDim = (XDimInfo*)call_data;
    vreal2 = &((XDim->Views)[XDim->oper2]);
    real = FI_GetLink("real", vreal2, XDim);
    imag = FI_GetLink("imag", vreal2, XDim);
    if(real < 0 || imag < 0)
      {
       strcpy(XDim->mess, "Link to REAL or IMAGinary part not found");
       WD_ViewportMessage(vreal2, XDim);
       return;
      }
    vreal2 = &((XDim->Views)[real]);
    if(vreal2->dataType)
      {
       WD_PopTypeCast(MY_DOUBLE, real, XDim);
       return;
      }
    vimag2 = &((XDim->Views)[imag]);
    if(vimag2->dataType)
      {
       WD_PopTypeCast(MY_DOUBLE, imag, XDim);
       return;
      }
    vreal1 = &((XDim->Views)[XDim->oper1]);
    real = FI_GetLink("real", vreal1, XDim);
    imag = FI_GetLink("imag", vreal1, XDim);
    vreal1 = &((XDim->Views)[real]);
    vimag1 = &((XDim->Views)[imag]);
    if(vreal1->actData == NULL || vimag1->actData == NULL)
      {
       strcpy(XDim->mess, "Couldn't process zero data field !");
       WD_ViewportMessage(vreal1, XDim);
       WD_ViewportMessage(vimag1, XDim);
       return;
      }
    if(vreal1->Width != vimag1->Width || vreal1->Height != vimag1->Height)
      {
       strcpy(XDim->mess, "Size of real and imag part must be the same !");
       WD_ViewportMessage(vreal1, XDim);
       WD_ViewportMessage(vimag1, XDim);
       return;
      }
    if(num)
      {
       ip = arg;
       x0 = -ip[0];
       y0 = -ip[1];
      }
    else
      {
       x0 = 0;
       y0 = 0;
      }
    vreal3 = WD_CopyViewport("(", "-", ")", vreal1, vreal2, XDim);
    vimag3 = WD_CopyViewport("(", "-", ")", vimag1, vimag2, XDim);
    if(vreal3 != NULL && vimag3 != NULL)
      {
       dr1 = vreal1->actData;
       di1 = vimag1->actData;
       dr3 = vreal3->actData;
       di3 = vimag3->actData;
       off = y0*vreal2->Width+x0;
       lr2 = vreal2->actData + off;
       li2 = vimag2->actData + off;
       y1 = y0;
       for(i=0; i<vreal1->Height; i++)
          {
           dr2 = lr2;
           di2 = li2;
           lr2 += vreal2->Width;
           li2 += vimag2->Width;
           x1 = x0;
           for(j=0; j < vreal1->Width; j++)
              {
               if(y1>=0 && y1 < vreal2->Height && x1 >=0 && x1<vreal2->Width)
                 {
                  *dr3++ = *dr1++ - *dr2++;
                  *di3++ = *di1++ - *di2++;
                 }
               else
                 {
                  *dr3++ = *dr1++;
                  *di3++ = *di1++;
        	  dr2++;
        	  di2++;
                 }
	       x1++;               
              }
           y1++;
          }
       CO_CalculateBases(vreal3);
       CO_CalculateBases(vimag3);
       FI_AddLink(vreal3->port, vimag3->port, "real", "imag", XDim);
       XT_Redraw(((XDim->myWidgets).Views[vreal3->port]).DArea,
       (XtPointer) vreal3, NULL);
       XT_Redraw(((XDim->myWidgets).Views[vimag3->port]).DArea,
       (XtPointer) vimag3, NULL);
      }
   }

/*
 * start complex sub
 */

void MA_StartCSub(Widget w, XtPointer client_data, XtPointer call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *vreal, *vimag;
    int real, imag;

    XDim = (XDimInfo*)client_data;
    vreal = &((XDim->Views)[XDim->port]);
    real = FI_GetLink("real", vreal, XDim);
    imag = FI_GetLink("imag", vreal, XDim);
    if(real < 0 || imag < 0)
      {
       strcpy(XDim->mess, "Link to REAL or IMAGinary part not found");
       WD_ViewportMessage(vreal, XDim);
       return;
      }
    vreal = &((XDim->Views)[real]);
    if(vreal->dataType)
      {
       WD_PopTypeCast(MY_DOUBLE, real, XDim);
       return;
      }
    vimag = &((XDim->Views)[imag]);
    if(vimag->dataType)
      {
       WD_PopTypeCast(MY_DOUBLE, imag, XDim);
       return;
      }
    if(vreal->actData == NULL || vimag->actData == NULL)
      {
       strcpy(XDim->mess, "Couldn't process zero data field !");
       WD_ViewportMessage(vreal, XDim);
       WD_ViewportMessage(vimag, XDim);
       return;
      }
    if(vreal->Width != vimag->Width || vreal->Height != vimag->Height)
      {
       strcpy(XDim->mess, "Size of real and imag part must be the same !");
       WD_ViewportMessage(vreal, XDim);
       WD_ViewportMessage(vimag, XDim);
       return;
      }
    WD_PopMath("Subtract from Complex Field", "Complex Field",
	real, MA_CSub, XDim);
   }

/*
 * multiply complex fields
 */

void MA_CMul(int num, void *arg, void *call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *vreal1, *vimag1, *vreal2, *vimag2, *vreal3, *vimag3;
    int i, j, *ip, x0, y0, x1, y1, off;
    int real, imag;
    double *dr1, *dr2, *dr3, *di1, *di2, *di3, re1, re2, im1, im2;
    double *li2, *lr2;

    XDim = (XDimInfo*)call_data;
    vreal2 = &((XDim->Views)[XDim->oper2]);
    real = FI_GetLink("real", vreal2, XDim);
    imag = FI_GetLink("imag", vreal2, XDim);
    if(real < 0 || imag < 0)
      {
       strcpy(XDim->mess, "Link to REAL or IMAGinary part not found");
       WD_ViewportMessage(vreal2, XDim);
       return;
      }
    vreal2 = &((XDim->Views)[real]);
    if(vreal2->dataType)
      {
       WD_PopTypeCast(MY_DOUBLE, real, XDim);
       return;
      }
    vimag2 = &((XDim->Views)[imag]);
    if(vimag2->dataType)
      {
       WD_PopTypeCast(MY_DOUBLE, imag, XDim);
       return;
      }
    vreal1 = &((XDim->Views)[XDim->oper1]);
    real = FI_GetLink("real", vreal1, XDim);
    imag = FI_GetLink("imag", vreal1, XDim);
    vreal1 = &((XDim->Views)[real]);
    vimag1 = &((XDim->Views)[imag]);
    if(vreal1->actData == NULL || vimag1->actData == NULL)
      {
       strcpy(XDim->mess, "Couldn't process zero data field !");
       WD_ViewportMessage(vreal1, XDim);
       WD_ViewportMessage(vimag1, XDim);
       return;
      }
    if(vreal1->Width != vimag1->Width || vreal1->Height != vimag1->Height)
      {
       strcpy(XDim->mess, "Size of real and imag part must be the same !");
       WD_ViewportMessage(vreal1, XDim);
       WD_ViewportMessage(vimag1, XDim);
       return;
      }
    if(num)
      {
       ip = arg;
       x0 = -ip[0];
       y0 = -ip[1];
      }
    else
      {
       x0 = 0;
       y0 = 0;
      }
    vreal3 = WD_CopyViewport("(", "*", ")", vreal1, vreal2, XDim);
    vimag3 = WD_CopyViewport("(", "*", ")", vimag1, vimag2, XDim);
    if(vreal3 != NULL && vimag3 != NULL)
      {
       dr1 = vreal1->actData;
       di1 = vimag1->actData;
       dr3 = vreal3->actData;
       di3 = vimag3->actData;
       off = y0*vreal2->Width+x0;
       lr2 = vreal2->actData + off;
       li2 = vimag2->actData + off;
       y1 = y0;
       for(i=0; i<vreal1->Height; i++)
          {
           dr2 = lr2;
           di2 = li2;
           lr2 += vreal2->Width;
           li2 += vimag2->Width;
           x1 = x0;
           for(j=0; j < vreal1->Width; j++)
              {
               if(y1>=0 && y1 < vreal2->Height && x1 >=0 && x1<vreal2->Width)
                 {
                  re1 = *dr1++;
                  im1 = *di1++;
                  re2 = *dr2++;
                  im2 = *di2++;
                  *dr3++ = re1*re2-im1*im2;
                  *di3++ = re1*im2+im1*re2;
                 }
               else
                 {
                  *dr3++ = *dr1++;
                  *di3++ = *di1++;
        	  dr2++;
        	  di2++;
                 }
	       x1++;               
              }
           y1++;
          }
       CO_CalculateBases(vreal3);
       CO_CalculateBases(vimag3);
       FI_AddLink(vreal3->port, vimag3->port, "real", "imag", XDim);
       XT_Redraw(((XDim->myWidgets).Views[vreal3->port]).DArea,
       (XtPointer) vreal3, NULL);
       XT_Redraw(((XDim->myWidgets).Views[vimag3->port]).DArea,
       (XtPointer) vimag3, NULL);
      }
   }

/*
 * start complex multiplication
 */

void MA_StartCMul(Widget w, XtPointer client_data, XtPointer call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *vreal, *vimag;
    int real, imag;

    XDim = (XDimInfo*)client_data;
    vreal = &((XDim->Views)[XDim->port]);
    real = FI_GetLink("real", vreal, XDim);
    imag = FI_GetLink("imag", vreal, XDim);
    if(real < 0 || imag < 0)
      {
       strcpy(XDim->mess, "Link to REAL or IMAGinary part not found");
       WD_ViewportMessage(vreal, XDim);
       return;
      }
    vreal = &((XDim->Views)[real]);
    if(vreal->dataType)
      {
       WD_PopTypeCast(MY_DOUBLE, real, XDim);
       return;
      }
    vimag = &((XDim->Views)[imag]);
    if(vimag->dataType)
      {
       WD_PopTypeCast(MY_DOUBLE, imag, XDim);
       return;
      }
    if(vreal->actData == NULL || vimag->actData == NULL)
      {
       strcpy(XDim->mess, "Couldn't process zero data field !");
       WD_ViewportMessage(vreal, XDim);
       WD_ViewportMessage(vimag, XDim);
       return;
      }
    if(vreal->Width != vimag->Width || vreal->Height != vimag->Height)
      {
       strcpy(XDim->mess, "Size of real and imag part must be the same !");
       WD_ViewportMessage(vreal, XDim);
       WD_ViewportMessage(vimag, XDim);
       return;
      }
    WD_PopMath("Multiply Complex Field", "by", real, MA_CMul, XDim);
   }

/*
 * divide complex fields
 */

void MA_CDiv(int num, void *arg, void *call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *vreal1, *vimag1, *vreal2, *vimag2, *vreal3, *vimag3;
    int i, j, *ip, x0, y0, x1, y1, off;
    int real, imag;
    double *dr1, *dr2, *dr3, *di1, *di2, *di3, re1, re2, im1, im2;
    double *li2, *lr2, val;
    Boolean Error = False;

    XDim = (XDimInfo*)call_data;
    vreal2 = &((XDim->Views)[XDim->oper2]);
    real = FI_GetLink("real", vreal2, XDim);
    imag = FI_GetLink("imag", vreal2, XDim);
    if(real < 0 || imag < 0)
      {
       strcpy(XDim->mess, "Link to REAL or IMAGinary part not found");
       WD_ViewportMessage(vreal2, XDim);
       return;
      }
    vreal2 = &((XDim->Views)[real]);
    if(vreal2->dataType)
      {
       WD_PopTypeCast(MY_DOUBLE, real, XDim);
       return;
      }
    vimag2 = &((XDim->Views)[imag]);
    if(vimag2->dataType)
      {
       WD_PopTypeCast(MY_DOUBLE, imag, XDim);
       return;
      }
    vreal1 = &((XDim->Views)[XDim->oper1]);
    real = FI_GetLink("real", vreal1, XDim);
    imag = FI_GetLink("imag", vreal1, XDim);
    vreal1 = &((XDim->Views)[real]);
    vimag1 = &((XDim->Views)[imag]);
    if(vreal1->actData == NULL || vimag1->actData == NULL)
      {
       strcpy(XDim->mess, "Couldn't process zero data field !");
       WD_ViewportMessage(vreal1, XDim);
       WD_ViewportMessage(vimag1, XDim);
       return;
      }
    if(vreal1->Width != vimag1->Width || vreal1->Height != vimag1->Height)
      {
       strcpy(XDim->mess, "Size of real and imag part must be the same !");
       WD_ViewportMessage(vreal1, XDim);
       WD_ViewportMessage(vimag1, XDim);
       return;
      }
    if(num)
      {
       ip = arg;
       x0 = -ip[0];
       y0 = -ip[1];
      }
    else
      {
       x0 = 0;
       y0 = 0;
      }
    off = y0*vreal2->Width+x0;
    lr2 = vreal2->actData + off;
    li2 = vimag2->actData + off;
    y1 = y0;
    for(i=0; i<vreal1->Height; i++)
       {
        dr2 = lr2;
        di2 = li2;
        lr2 += vreal2->Width;
        li2 += vimag2->Width;
        x1 = x0;
        for(j=0; j < vreal1->Width; j++)
           {
            if(y1>=0 && y1 < vreal2->Height && x1 >=0 && x1<vreal2->Width)
              {
               re2 = *dr2;
               im2 = *di2;
               if(re2*re2+im2*im2 == 0)
                 {
                  i = vreal1->Height;
                  Error = True;
                  break;
                 }
              }
            dr2++;
            di2++;
	    x1++;               
           }
        y1++;
       }
    if(Error)
      {
       strcpy(XDim->error,
	     "To divide all selected field elements must be != 0");
       WD_XDimError(XDim, NULL);
       vreal3 = NULL;
       vimag3 = NULL;
      }
    else
      {
       vreal3 = WD_CopyViewport("(", "/", ")", vreal1, vreal2, XDim);
       vimag3 = WD_CopyViewport("(", "/", ")", vimag1, vimag2, XDim);
      }
    if(vreal3 != NULL && vimag3 != NULL)
      {
       dr1 = vreal1->actData;
       di1 = vimag1->actData;
       dr3 = vreal3->actData;
       di3 = vimag3->actData;
       lr2 = vreal2->actData + off;
       li2 = vimag2->actData + off;
       y1 = y0;
       for(i=0; i<vreal1->Height; i++)
          {
           dr2 = lr2;
           di2 = li2;
           lr2 += vreal2->Width;
           li2 += vimag2->Width;
           x1 = x0;
           for(j=0; j < vreal1->Width; j++)
              {
               if(y1>=0 && y1 < vreal2->Height && x1 >=0 && x1<vreal2->Width)
                 {
                  re1 = *dr1++;
                  im1 = *di1++;
                  re2 = *dr2++;
                  im2 = *di2++;
                  val = 1.0/(re2*re2+im2*im2);
                  *dr3++ = (re1*re2+im1*im2)*val;
                  *di3++ = (im1*re2-re1*im2)*val;
                 }
               else
                 {
                  *dr3++ = *dr1++;
                  *di3++ = *di1++;
                 }
               dr2++;
               di2++;
	       x1++;               
              }
           y1++;
          }
       CO_CalculateBases(vreal3);
       CO_CalculateBases(vimag3);
       FI_AddLink(vreal3->port, vimag3->port, "real", "imag", XDim);
       XT_Redraw(((XDim->myWidgets).Views[vreal3->port]).DArea,
       (XtPointer) vreal3, NULL);
       XT_Redraw(((XDim->myWidgets).Views[vimag3->port]).DArea,
       (XtPointer) vimag3, NULL);
      }
   }

/*
 * start complex division
 */

void MA_StartCDiv(Widget w, XtPointer client_data, XtPointer call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *vreal, *vimag;
    int real, imag;

    XDim = (XDimInfo*)client_data;
    vreal = &((XDim->Views)[XDim->port]);
    real = FI_GetLink("real", vreal, XDim);
    imag = FI_GetLink("imag", vreal, XDim);
    if(real < 0 || imag < 0)
      {
       strcpy(XDim->mess, "Link to REAL or IMAGinary part not found");
       WD_ViewportMessage(vreal, XDim);
       return;
      }
    vreal = &((XDim->Views)[real]);
    if(vreal->dataType)
      {
       WD_PopTypeCast(MY_DOUBLE, real, XDim);
       return;
      }
    vimag = &((XDim->Views)[imag]);
    if(vimag->dataType)
      {
       WD_PopTypeCast(MY_DOUBLE, imag, XDim);
       return;
      }
    if(vreal->actData == NULL || vimag->actData == NULL)
      {
       strcpy(XDim->mess, "Couldn't process zero data field !");
       WD_ViewportMessage(vreal, XDim);
       WD_ViewportMessage(vimag, XDim);
       return;
      }
    if(vreal->Width != vimag->Width || vreal->Height != vimag->Height)
      {
       strcpy(XDim->mess, "Size of real and imag part must be the same !");
       WD_ViewportMessage(vreal, XDim);
       WD_ViewportMessage(vimag, XDim);
       return;
      }
    WD_PopMath("Divide Complex Field", "by", real, MA_CMul, XDim);
   }

/*
 * Calculate FFT, l2 <= log2(len)+1
 * (Field length = 2^l2, fact is used to normalize, typically
 * fact = 1/N or fact = 1/sqr(N) etc.)
 * if fact <0 inv. FFT assumed
 * step gives distance between two values (p. ex. xy fft)
 */

void MA_CalcFFT(double *real, double *imag, int l2, int step, double fact)
   {
    int len, len2, le, le1, i, j, l, bitRev, m1, m2, sm2, inv, s1, s2;
    double a, wr, wi, ur, ui, tr, ti;

    if(real != NULL && imag !=NULL)
      {
       len = 1<<l2;
       len2 = len >>1;
       sm2 = 1<<(l2-1);
       s1 = 0;
       for(i=0; i <len; i++)
          {
           bitRev = 0;
	   m1 = 1;
	   m2 = sm2;
           while(m2>=m1)
              {
               if(i&m1)
                 bitRev |= m2;
               if(i&m2)
                 bitRev |= m1;
               m2 >>= 1;
               m1 <<= 1;
              }
           if(bitRev > i)
             {
              s2 = bitRev*step;
              a = real[s2];
	      real[s2] = real[s1];
	      real[s1] = a;
	      a = imag[s2];
	      imag[s2] = imag[s1];
	      imag[s1] = a;
	     }
           s1 += step;
          }
       if(fact < 0)
         {
          inv = False;
          fact = -fact;
         }
       else
         inv = True;
       le = 1;
       for(l=0; l<l2; l++)
         {
          le1 = le;
          le <<= 1;
          a = PI/le1;
          wr = cos(a);
          wi = sin(a);
          if(inv)
            wi = -wi;
          m1 = le1*step;
          ur = 1.0;
          ui = 0.0;
          for(j=0; j<le1; j++)
             {
              for(i=j; i<len; i+=le)
                 {
                  s1 = i*step;
                  s2 = (i+le1)*step;
                  tr = real[s2]*ur-imag[s2]*ui;
                  ti = real[s2]*ui+imag[s2]*ur;
                  real[s2] = real[s1]-tr;
                  imag[s2] = imag[s1]-ti;
                  real[s1] += tr;
                  imag[s1] += ti;
                 }
              a = ur*wr-ui*wi;
              ui = ur*wi+ui*wr;
              ur = a;
             }
         }
       s1 = 0;
       for(i=0; i< len; i++)
         {
          real[s1] *= fact;
          imag[s1] *= fact;
          s1 += step;
         }
      }
   }

/*
 * start FFT
 */

void MA_StartFFT(Widget w, XtPointer client_data, XtPointer call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *vreal, *vimag, *vimag2, *vreal2;
    int real, imag, newWidth, newHeight, l2w, l2h, port, i, j, fft;
    double *dr, *di, factX, factY, x, y;

    XDim = (XDimInfo*)client_data;
    vreal = &((XDim->Views)[XDim->port]);
    real = FI_GetLink("real", vreal, XDim);
    imag = FI_GetLink("imag", vreal, XDim);
    if(real < 0 || imag < 0)
      {
       if(vreal->dataType)
	 {
	  WD_PopTypeCast(MY_DOUBLE, XDim->port, XDim);
	  return;
	 }
       strcpy(XDim->mess,
	"Link to REAL or IMAGinary part not found => used as real");
       WD_ViewportMessage(vreal, XDim);
       vimag = NULL;
      }
    else
      {
       vreal = &((XDim->Views)[real]);
       if(vreal->dataType)
	 {
	  WD_PopTypeCast(MY_DOUBLE, real, XDim);
	  return;
	 }
       vimag = &((XDim->Views)[imag]);
       if(vimag->dataType)
	 {
	  WD_PopTypeCast(MY_DOUBLE, imag, XDim);
	  return;
	 }
      }
    if(vreal->actData == NULL || (vimag != NULL && vimag->actData == NULL))
      {
       strcpy(XDim->mess, "Couldn't process zero data field !");
       WD_ViewportMessage(vreal, XDim);
       if(vimag != NULL)
         WD_ViewportMessage(vimag, XDim);
       return;
      }
    if(vimag != NULL &&
	(vreal->Width != vimag->Width || vreal->Height != vimag->Height))
      {
       strcpy(XDim->mess, "Size of real and imag part must be the same !");
       WD_ViewportMessage(vreal, XDim);
       WD_ViewportMessage(vimag, XDim);
       return;
      }
    newWidth = 1;
    newHeight = 1;
    l2w = 0;
    l2h = 0;
    while(newWidth < vreal->Width || newHeight < vreal->Height)
       {
        if(newWidth < vreal->Width)
          {
           newWidth <<=1;
           l2w++;
          }
        if(newHeight < vreal->Height)
          {
           newHeight <<=1;
           l2h++;
          }
       }
    port = XDim->port;
    WD_CreateViewport((XDim->myWidgets).NewView, XDim, NULL);
    if(port != XDim->port)
      {
       if(CO_NewData(newWidth, newHeight, MY_DOUBLE, XDim))
	 {
	  WD_XDimError(XDim, NULL);
	  XT_DestroyViewport(NULL,&(XDim->port), NULL);
	  return;
	 }
      }
    else
      return;
    vreal2 = &((XDim->Views)[XDim->port]);
    WD_SetTitle(XDim->port, vreal, "real(FFT(", "))", XDim);
    port = XDim->port;
    WD_CreateViewport((XDim->myWidgets).NewView, XDim, NULL);
    if(port != XDim->port)
      {
       if(CO_NewData(newWidth, newHeight, MY_DOUBLE, XDim))
	 {
	  WD_XDimError(XDim, NULL);
	  XT_DestroyViewport(NULL,&(XDim->port), NULL);
	  XT_DestroyViewport(NULL,&(vreal2->port), NULL);
	  return;
	 }
      }
    else
      {
       XT_DestroyViewport(NULL,&(vreal2->port), NULL);
       return;
      }
    vimag2 = &((XDim->Views)[XDim->port]);
    WD_SetTitle(XDim->port, vreal, "imag(FFT(", "))", XDim);
    if(newWidth > 1)
      factX = ((double)(vreal->Width-1))/(newWidth-1);
    else
      factX = 1.0;
    if(newHeight > 1)
      factY = ((double)(vreal->Height-1))/(newHeight-1);
    else
      factY = 1.0;
    vreal2->FactX = factX*vreal->FactX;
    vreal2->FactY = factY*vreal->FactY;
    vimag2->FactX = factX*vreal->FactX;
    vimag2->FactY = factY*vreal->FactY;
    /*
     * interpolate data, to use spline interpolation
     * use spline from the data menu and set it to a size=2^len
     * the minus sign is used to shift the data to the middle
     */
    if(w==(XDim->myWidgets).ViewMenus.fft)
      fft = True;
    else
      fft = False;
    dr = vreal2->actData;
    di = vimag2->actData;
    y = 0.0;
    if(vimag != NULL)
      {
       for(i=0; i<newHeight; i++)
	  {
	   x = 0.0;
	   for(j=0; j<newWidth; j++)
	      {
	       if(fft && ((i+j)&1))
	         {
	          *dr++ = -CO_Inter(x,y,vreal);
	          *di++ = -CO_Inter(x,y,vimag);
	         }
	       else
	         {
	          *dr++ = CO_Inter(x,y,vreal);
	          *di++ = CO_Inter(x,y,vimag);
	         }
	       x += factX;
	      }
	   y += factY;
	  }
       }
     else
      {
       for(i=0; i<newHeight; i++)
	  {
	   x = 0.0;
	   for(j=0; j<newWidth; j++)
	      {
	       if(fft && ((i+j)&1))
	         {
		  *dr++ = -CO_Inter(x,y,vreal);
		  *di++ = 0.0;
		 }
	       else
	         {
		  *dr++ = CO_Inter(x,y,vreal);
		  *di++ = 0.0;
		 }
	       x += factX;
	      }
	   y += factY;
	  }
       }
    /*
     * define your normalization factors here
     */
    if(fft)
      {
       factX = 1.0/newWidth;
       factY = 1.0/newHeight;
      }
    else
      {
      /* inverse FFT */
       factX = -1.0;
       factY = -1.0;
      }
    dr = vreal2->actData;
    di = vimag2->actData;
    for(i=0; i<newHeight; i++)
       {
	MA_CalcFFT(dr, di, l2w, 1, factX);
	dr += newWidth;
	di += newWidth;
       }
    dr = vreal2->actData;
    di = vimag2->actData;
    for(i=0; i<newWidth; i++)
       {
	MA_CalcFFT(dr, di, l2h, newWidth, factY);
	dr++;
	di++;
       }
    if(!fft)
      {
       dr = vreal2->actData;
       di = vimag2->actData;
       for(i=0; i<newHeight;i++)
          {
           for(j=0; j<newWidth;j++)
              {
               if((i+j)&1)
                 {
                  *dr = -*dr;
                  *di = -*di;
                 }
               dr++;
               di++;
              }
          }
      }
    CO_CalculateBases(vreal2);
    CO_CalculateBases(vimag2);
    FI_AddLink(vreal2->port, vimag2->port, "real", "imag", XDim);
    XT_Redraw(((XDim->myWidgets).Views[vreal2->port]).DArea,
    (XtPointer) vreal2, NULL);
    XT_Redraw(((XDim->myWidgets).Views[vimag2->port]).DArea,
    (XtPointer) vimag2, NULL);
   }

/*
 * Create linear Low/Highpass filters
 */


void MA_HiLoPass(int num, void *arg, void *call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *view, *vreal, *vimag;
    int real, imag, cut, ord, v0, v1, i, j, delx, dely;
    double *dr, *di, re, D0, fact;
    myString *vals;
    myString name;

    XDim = (XDimInfo*)call_data;
    vals = (myString *) arg;
    if(!sscanf(vals[0], "%le", &D0) || D0 <= 0.0)
      {
       strcpy(XDim->error, "Wrong Filter Radius!!");
       WD_XDimError(XDim, NULL);
       return;
      }
    if(num == 2)
      {
       ord = 1;
       if(!sscanf(vals[1], "%d", &ord) || ord <=0)
	 {
	  strcpy(XDim->error, "Error in Order Number!!");
	  WD_XDimError(XDim, NULL);
	  return;
	 }
      }
    if(XDim->port >=0)
      {
       view = &((XDim->Views)[XDim->port]);
       real = FI_GetLink("real", view, XDim);
       imag = FI_GetLink("imag", view, XDim);
       if(real < 0 || imag < 0)
         {
          strcpy(XDim->mess, "Link to REAL or IMAGinary part not found");
          WD_ViewportMessage(view, XDim);
         }
       else
         {
          view = &((XDim->Views)[real]);
	 }          
       if(view->Width <= 0 || view->Height <= 0)
         {
          strcpy(XDim->mess, "Couldn't process zero data field !");
          WD_ViewportMessage(view, XDim);
          return;
         }
       switch (XDim->aView.filter)
         {
          default:
          case IdealL:
          case IdealH:
	       strcpy(name, "Ideal Filter for");
	       break;
	  case ButwL:
	  case ButwH:
	       strcpy(name, "Butterworth Filter for");
	       break;
	  case GaussL:
	  case GaussH:
	       strcpy(name, "Gauss Filter for");
	 }
       vreal = WD_CopyViewport(name, "", " (real)", view, NULL, XDim);
       vimag = WD_CopyViewport(name, "", " (imag)", view, NULL, XDim);
       if(vreal != NULL && vimag != NULL)
         {
          real = vreal->port;
          imag = vimag->port;
	  dr = vreal->actData;
	  di = vimag->actData;
	  D0 *= view->Width;
	  D0 /= 100.0;
	  dely = view->Height >> 1;
	  delx = view->Width >> 1;
	  switch(XDim->aView.filter)
	     {
	      default:
	      case IdealL:
	        cut = (int) (D0*D0);
		for(i = 0; i<view->Height; i++)
        	   {
        	    v0 = (i-dely);
        	    v0 *= v0;
        	    for(j=0; j<view->Width; j++)
        	       {
			*di++ = 0.0;
			v1 = j-delx;
			v1 *= v1;
			if(v0+v1<cut)
			  *dr++ = 1.0;
			else
			  *dr++ = 0.0;
		       }
        	   }
        	break;
              case IdealH:
	        cut = (int) (D0*D0);
		for(i = 0; i<view->Height; i++)
        	   {
        	    v0 = (i-dely);
        	    v0 *= v0;
        	    for(j=0; j<view->Width; j++)
        	       {
			*di++ = 0.0;
			v1 = j-delx;
			v1 *= v1;
			if(v0+v1<cut)
			  *dr++ = 0.0;
			else
			  *dr++ = 1.0;
		       }
        	   }
        	break;
              case GaussL:
	        fact =  log(0.5)/(D0*D0);
		for(i = 0; i<view->Height; i++)
        	   {
        	    v0 = (i-dely);
        	    v0 *= v0;
        	    for(j=0; j<view->Width; j++)
        	       {
			*di++ = 0.0;
			v1 = j-delx;
			v1 *= v1;
			re = exp(fact*((double)(v0+v1)));
			*dr++ = re;
		       }
        	   }
        	break;
              case GaussH:
	        fact =  log(0.5)/(D0*D0);
		for(i = 0; i<view->Height; i++)
        	   {
        	    v0 = (i-dely);
        	    v0 *= v0;
        	    for(j=0; j<view->Width; j++)
        	       {
			*di++ = 0.0;
			v1 = j-delx;
			v1 *= v1;
			re = 1-exp(fact*((double)(v0+v1)));
			*dr++ = re;
		       }
        	   }
        	break;
              case ButwL:
	        fact =  1.0/(D0*D0);
		for(i = 0; i<view->Height; i++)
        	   {
        	    v0 = (i-dely);
        	    v0 *= v0;
        	    for(j=0; j<view->Width; j++)
        	       {
			*di++ = 0.0;
			v1 = j-delx;
			v1 *= v1;
			re = pow(fact*((double)(v0+v1)), (double)ord);
			*dr++ = 1.0/(1.0+re);
		       }
        	   }
        	break;
              case ButwH:
	        fact =  D0*D0;
		for(i = 0; i<view->Height; i++)
        	   {
        	    v0 = (i-dely);
        	    v0 *= v0;
        	    for(j=0; j<view->Width; j++)
        	       {
			*di++ = 0.0;
			v1 = j-delx;
			v1 *= v1;
			if(v0+v1 == 0)
			  re = 0.0;
			else
			  re = 1.0/(1.0+pow(fact/((double)(v0+v1)),
				(double)ord));
			*dr++ = re;
		       }
        	   }
        	break;
             }
          CO_CalculateBases(vreal);
          CO_CalculateBases(vimag);
          FI_AddLink(real, imag, "real", "imag", XDim);
	  XT_Redraw(((XDim->myWidgets).Views[real]).DArea,
	  (XtPointer) vreal, NULL);
	  XT_Redraw(((XDim->myWidgets).Views[imag]).DArea,
	  (XtPointer) vimag, NULL);
	 }
      }
   }

/*
 * get radius of Hi-/Lo Filter
 */

void MA_StartHiLo(Widget w, XtPointer client_data, XtPointer call_data)
   {
    extern XDimInfo XDim;

    XDim.aView.filter = (enum XDFilter) client_data;
    if(XDim.aView.filter == ButwL || XDim.aView.filter == ButwH)
       WD_PopInput2("Filter Radius: (in % from Width)", "50", 5, 10,
	"Order:", "1", 5, 5, MA_HiLoPass, &XDim);
    else
       WD_PopInput("Filter Radius (in % from Width):", "50", 5, 10,
	MA_HiLoPass, &XDim);
   }

/*
 * Sort program for qsort
 */

int MA_SortVal(const void *v1, const void *v2)
   {
    int ret;
    const double *val1=v1, *val2=v2;
    double val;
    
    val = *val1-*val2;
    if(val==0.0)
      ret =0;
    else
      ret = (val>0.0 ? 1 : -1);
    return(ret);
   }

/*
 * Median Filter
 */

void MA_Median(int num, void *arg, void *call_data)
   {
    XDimInfo *XDim;
    ViewportInfo *view1, *view2;
    myString *vals;
    double *dp1, *dp2, *vec;
    int Width, Height, x, y, i, j, k, dely, vx, vy, size;

    XDim = (XDimInfo*)call_data;
    if(XDim->port >=0)
      {
       view1 = &((XDim->Views)[XDim->port]);
       if(view1->dataType)
	  WD_PopTypeCast(MY_DOUBLE, XDim->port, XDim);
       if(view1->actData == NULL)
         return;
       vals = (myString *) arg;
       if(!sscanf(vals[0], "%d", &Width) || Width <= 0)
	 {
	  strcpy(XDim->error, "Error in Filter Width!!");
	  WD_XDimError(XDim, NULL);
	  return;
	 }
       if(!sscanf(vals[1], "%d", &Height) || Height <=0)
	 {
	  strcpy(XDim->error, "Error in Filter Height!!");
	  WD_XDimError(XDim, NULL);
	  return;
	 }
       /*
        * Width and Height must be uneven
        */
       if(Width &1 == 0 || Height &1 == 0)
         {
          strcpy(XDim->error, "Width and Height are supposed to be uneven!!");
          WD_ViewportMessage(view1, XDim);
          Width |= 1;
          Height |= 1;
         }
       size = Width*Height;
       /*
	* create sorting vector
	*/
       vec = (double*)malloc(size*sizeof(double));
       if(vec == NULL)
         {
          sprintf(XDim->error, "Couldn't allocate %d number of bytes !!!",
		size*(int)sizeof(double));
	  WD_XDimError(XDim, NULL);
	  return;
         }
       view2 = WD_CopyViewport("Median(", "", ")", view1, NULL, XDim);
       if(view2)
         {
	  dp1 = view1->actData;
	  dp2 = view2->actData;
	  /*
	   * Width and Height uneven => symmetric to zero
	   */
	  Width >>= 1;
	  Height >>= 1;
	  for(y=0; y<view1->Height; y++)
	     {
	      for(x=0; x<view1->Width; x++)
	         {
	          k = 0;
	          for(i= -Height; i<=Height; i++)
	             {
	              vy= y+i;
	              dely = i*view1->Width;
	              for(j= -Width; j<=Width; j++)
	                 {
	                  vx =x+j;
	                  if(vy <0 || vy >= view1->Height||
				vx <0 || vx >= view1->Width)
	                    vec[k++] = *dp1;
	                  else
	                    vec[k++] = *(dp1+dely+j);
	                 }
	             }
	          dp1++;
	          qsort(vec, size, sizeof(double), MA_SortVal);
	          *dp2++ = vec[(size>>1)];
	         }
	     }
          CO_CalculateBases(view2);
	  XT_Redraw(((XDim->myWidgets).Views[XDim->port]).DArea,
	  (XtPointer) view2, NULL);
	 }
      }
   }

/*
 * Matrix Filter
 */

void MA_Matrix(int num, void *arg, void *call_data)
   {
    XDimInfo *XDim;
    XDimEditMatrix *XEM;
    ViewportInfo *view1, *view2;
    double *dp1, *dp2, *dm, Sum, Fact, div;
    int Width, Height, x, y, i, j, dely, vx, vy;

    XDim = (XDimInfo*)call_data;
    if(XDim->port >=0)
      {
       view1 = &((XDim->Views)[XDim->port]);
       if(view1->dataType)
	  WD_PopTypeCast(MY_DOUBLE, XDim->port, XDim);
       if(view1->actData == NULL)
         return;
       view2 = WD_CopyViewport("Matrix(", "", ")", view1, NULL, XDim);
       if(view2)
         {
          XEM = &((XDim->myWidgets).EditMatrix);
	  dp1 = view1->actData;
	  dp2 = view2->actData;
	  /*
	   * Width and Height uneven => symmetric to zero
	   */
	  Width = (XEM->Width)>>1;
	  Height = (XEM->Height)>>1;
          if(XEM->div != 0.0)
            div = 1.0/XEM->div;
          else
            div = 1.0;
	  /*
	   * if outside use factor of middle value
	   */
	  Fact = (XEM->Values)[Width*(XEM->Width)+Height];
	  for(y=0; y<view1->Height; y++)
	     {
	      for(x=0; x<view1->Width; x++)
	         {
	          dm = XEM->Values;
	          Sum = 0.0;
	          for(i= -Height; i<=Height; i++)
	             {
	              vy= y+i;
	              dely = i*view1->Width;
	              for(j= -Width; j<=Width; j++)
	                 {
	                  vx =x+j;
	                  if(vy <0 || vy >= view1->Height||
				vx <0 || vx >= view1->Width)
	                    Sum += Fact * (*dp1);
	                  else
	                    Sum +=  (*dm) * *(dp1+dely+j);
	                  dm++;
	                 }
	             }
	          *dp2++ = (div*Sum)+ XEM->bias;
	          dp1++;
	         }
	     }
          CO_CalculateBases(view2);
	  XT_Redraw(((XDim->myWidgets).Views[XDim->port]).DArea,
	  (XtPointer) view2, NULL);
	 }
      }
   }


/*
 * Pop up Matrix
 */

void MA_PopMatrix(int num, void *arg, void *call_data)
   {
    XDimInfo *XDim;
    XDimEditMatrix *XEM;
    ViewportInfo *view1;
    myString *vals;
    char val[32];
    int Width, Height;

    XDim = (XDimInfo*)call_data;
    if(XDim->port >=0)
      {
       XEM = &((XDim->myWidgets).EditMatrix);
       view1 = &((XDim->Views)[XDim->port]);
       if(view1->actData == NULL)
         return;
       vals = (myString *) arg;
       if(!sscanf(vals[0], "%d", &Width) || Width <= 0)
	 {
	  strcpy(XDim->error, "Error in Filter Width!!");
	  WD_XDimError(XDim, NULL);
	  return;
	 }
       if(!sscanf(vals[1], "%d", &Height) || Height <=0)
	 {
	  strcpy(XDim->error, "Error in Filter Height!!");
	  WD_XDimError(XDim, NULL);
	  return;
	 }
       /*
        * Width and Height must be uneven
        */
       if((Width &1) == 0 || (Height &1) == 0)
         {
          strcpy(XDim->mess, "Width and Height are supposed to be uneven!!");
          WD_ViewportMessage(view1, XDim);
          Width |= 1;
          Height |= 1;
         }
       sprintf(val, "%d", Width*Height);
       if(WD_CreateEditMatrix("1",val,"0", Width, Height, -(Width>>1),
		-(Height>>1), XDim))
         {
	  XtManageChild(XEM->Shell);
	 }
       else
         {
	  strcpy(XDim->error, "Failed to allocate Edit Matrix!!");
	  WD_XDimError(XDim, NULL);
	  return;
         }
      }
   }

/*
 * get width and height of neightbourhood
 */

void MA_StartOtherF(Widget w, XtPointer client_data, XtPointer call_data)
   {
    XDimInfo *XDim;
    
    XDim = (XDimInfo*) client_data;
    if(w==(XDim->myWidgets).ViewMenus.median)
       WD_PopInput2("Median Filter Width (uneven):", "3", 4, 4,
	"Height (uneven):", "3", 4, 4, MA_Median, XDim);
    else
      {
       (XDim->myWidgets).EditMatrix.RetFunc = MA_Matrix;
       WD_PopInput2("Matrix Filter Width (uneven):", "3", 4, 4,
	"Height (uneven):", "3", 4, 4, MA_PopMatrix, XDim);
      }
   }
