/* some of Mats' routines written by W.W. , adapted for ATC version 12/6/96 */
 /* 11/20/97 - some mods in ufourn, usvbksb, usvdcmp to make them safer,
 they were assuming F*4 arrays */
 /*call user defined routines from here            W.W.......1992-11-19*/
 
 /* To check variable class & type:
		 sym[?].class:	8-scalar ptr;	2-string
				 1-scalar;	4 array
		 sym[?].type:	0-byte	1-int	2-long
 
 				3-float		4-double  */
 
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include "ana_structures.h"
#if defined(SOLARIS)
#include <unistd.h>
#else
#include <sys/unistd.h>
#endif
 extern	struct sym_desc sym[];
 extern	struct sym_list		*subr_sym_list[];
 extern	int	ana_type_size[];
 extern	int	vfix_top, num_ana_subr, next_user_subr_num;
 extern	int	lastmin_sym, lastmax_sym;
 struct  ahead   *h;
 union	types_ptr { byte *b; short *w; int *l; float *f; double *d;};		
 union	types_ptr q1,q2,q3,q4,q5;
 int ni_save=0, nj_save=0;
 /* note - a lot of stuff that was commented out was removed here, see Mats'
 original if of interest */
 /*------------------------------------------------------------------------- */
 /* pgm format read. pgm format is used when storing videk images(videk.c).
		 fp=fopen("q0.pgm","w");
	sprintf(head, "P5\n%d %d\n%d\n", w, h, 255);
	 fwrite(head, 1, strlen(head), fp);
	 fwrite( ((char *)buf_virt_addr[0]), 1, w*h, fp);
 --Images like--------------------
 P5
 256 256
 255
 PKGJDC:?A;;:;54/.-+'(&'%#%'$!'%$$! $#   <--binary(8 bit) data
 ----------------------
*/
int ana_pgmread(narg,ps)
 /* PROCEDURE -- Called by: PGMREAD,VAR,'FILENAME'*/
 int     narg, ps[];
 {
 int    dim[2],iq, n, nb,i,pixelsizecode;
 struct ahead   *h;
 union  types_ptr q1;
 char buf[50],*name;
 FILE *fopen(),*fin;
 char    *ulib_path;

         /* first arg is the variable to load, second is name of file */
 if ( sym[ ps[1] ].class != 2 ) return execute_error(70);

 /* try to open the file */
  name = (char *) sym[ps[1] ].spec.array.ptr;
 if ((fin=fopen(name,"r")) == NULL)
 {
 sprintf(buf,"%s%s",ulib_path,name);
  if ((fin=fopen(buf,"r")) == NULL) return execute_error(84);
 }

        fscanf(fin,"%s %d %d %d",buf,&dim[0],&dim[1],&i);
 /*      printf("\n%s %d %d %d\n",buf,dim[0],dim[1],i); */

        switch(i){
                case 0xff:      { pixelsizecode=0; break; }  /* 1 byte */
                case 0xffff:    { pixelsizecode=1; break; }  /* 2 bytes */
                case 0xffffffff:{ pixelsizecode=2; break;}   /* 4 bytes */
                }
 /* compute size of array in bytes*/
 nb=dim[0]*dim[1]*ana_type_size[pixelsizecode];

         /* create the output array(data length: 1 byte & 2 dimensions) */
 iq = ps[0];
 if ( redef_array( iq, pixelsizecode, 2, dim) != 1 ) {
         fclose(fin);
         return -1; }

 h = (struct ahead *) sym[iq].spec.array.ptr;
 q1.l = (int *) ((char *)h + sizeof(struct ahead));
                                                 /* just read into the array */
 fseek(fin,1,SEEK_CUR);          /* move forward to overcome '0a' (==LF) byte */
 if (fread(q1.l, 1, nb, fin) != nb) {
         execute_error(92);  fclose(fin);  return -1; }

 /* common section again */
 fclose(fin);
 return 1;
 }
 /*------------------------------------------------------------------------- */
	/* call fourn.f. fourn(float data[], int nn[], int ndim, int sign) */
int ufourn(narg,ps)
 int narg, ps[];
 {
 /*	extern fourn(); */
	 int iq, i1, i2, data_type;
 
 if(narg < 3) {printf("\nNot enough parameters to call fourn()\n");
	 return -1;              }
 
  data_type = sym[ps[0]].type;
  if(data_type != 3) {
	  printf("data_type....%d\n",data_type);
	  return -1;               /* not float type */
  }
	 /* array ptr of float data[] --> q1.f */
 iq=ps[0];
 h = (struct ahead *) sym[iq].spec.array.ptr;
 q1.f = (float *) ((char *) h + sizeof( struct ahead));
 
	 /* array ptr of int nn[] --> p2.l */
 iq=ps[1];
 h = (struct ahead *) sym[iq].spec.array.ptr;
 q2.l = (int *) ((char *) h + sizeof( struct ahead));
 
	 /* value of ndim --> i1 */
 if (int_arg_stat(ps[2], &i1) != 1) return -1;
 
	 /* value of sign --> i2 */
 if (int_arg_stat(ps[3], &i2) != 1) return -1;
 
	 /* call fortran routine fourn.f */
	 fourn_(q1.f, q2.l, &i1, &i2);
 /*	fourn(q1.f-1, q2.l-1, i1, i2);*/	/* call fourn.c, C routime */
 
 eend: return 1;
 }
 /*------------------------------------------------------------------------- */
	/* call svbksb.f-- svbksb(u, w, v, m, n, mp, np, b, x)
		float u[m][n];	float w[n];	float v[n][n];
		int m, n, mp, np;
		float b[m];	float x[n];
	*/
int usvbksb(narg,ps)
 int narg, ps[];
 {
	 extern SVBKSB();
	 int i1,i2,i3,i4,iq;
 
 if(narg < 8) {printf("\nNot enought parameters to call usvbksb()\n");
		 return -1;		}
 
	 /* float u[][] ptr --> q1.f */
 iq = ana_float(1, &ps[0]);
 h = (struct ahead *) sym[iq].spec.array.ptr;
 q1.f = (float *) ((char *) h + sizeof( struct ahead));
 
	 /* float w[] ptr --> q2.f */
 iq = ana_float(1, &ps[1]);
 h = (struct ahead *) sym[iq].spec.array.ptr;
 q2.f = (float *) ((char *) h + sizeof( struct ahead));
 
	 /* float v[][]ptr --> q3.f */
 iq = ana_float(1, &ps[2]);
 h = (struct ahead *) sym[iq].spec.array.ptr;
 q3.f = (float *) ((char *) h + sizeof( struct ahead));
 
	 /* i1 == int m */
 if (int_arg_stat(ps[3], &i1) != 1) return -1;
 
	 /* i2 == int n */
 if (int_arg_stat(ps[4], &i2) != 1) return -1;
 
	 /* i3 == int mp */
 if (int_arg_stat(ps[5], &i3) != 1) return -1;
 
	 /* i4 == int np */
 if (int_arg_stat(ps[6], &i4) != 1) return -1;
 
	 /* float b[] ptr --> q4.f */
 iq = ana_float(1, &ps[7]);
 h = (struct ahead *) sym[iq].spec.array.ptr;
 q4.f = (float *) ((char *) h + sizeof( struct ahead));
       
	 /* float x[] ptr --> q5.f */
 iq = ana_float(1, &ps[8]);
 h = (struct ahead *) sym[iq].spec.array.ptr;
 q5.f = (float *) ((char *) h + sizeof( struct ahead));
 
	 /* fortran subroutine SVBKSB() call */
 svbksb_(q1.f, q2.f, q3.f, &i1, &i2, &i3, &i4, q4.f, q5.f);
 
 return 1;
 }
 /*------------------------------------------------------------------------- */
	/* call svdcmp.f -- svdcmp(a, m, n, mp, np, w, v)
		float a[][], w[], v[][];
		int m, n, mp, np;
	*/
int usvdcmp(narg,ps)
 int narg, ps[];
 {
	 extern svdcmp();
	 int i1,i2,i3,i4,iq;
 
 if(narg < 6) {printf("\nNo enought parameters to call svdcmp()\n");
	 return -1;              }
 
	 /* get ptr of a[][] -- q1.f */
 iq = ana_float(1, &ps[0]);
 h = (struct ahead *) sym[iq].spec.array.ptr;
 q1.f = (float *) ((char *) h + sizeof( struct ahead));
 
	 /* i1 == int m */
 if (int_arg_stat(ps[1], &i1) != 1) return -1;
 
	 /* i2 == int n */
 if (int_arg_stat(ps[2], &i2) != 1) return -1;
			    
	 /* i3 == mp */
 if (int_arg_stat(ps[3], &i3) != 1) return -1;
 
	 /* i4 == int np */
 if (int_arg_stat(ps[4], &i4) != 1) return -1;
 
	 /* get ptr of w[] -- q2.f */
 iq = ana_float(1, &ps[5]);
 h = (struct ahead *) sym[iq].spec.array.ptr;
 q2.f = (float *) ((char *) h + sizeof( struct ahead));
     
	 /* get ptr of v[][] -- q3.f */
 iq = ana_float(1, &ps[6]);
 h = (struct ahead *) sym[iq].spec.array.ptr;
 q3.f = (float *) ((char *) h + sizeof( struct ahead));
 
	 /* call SVDCMP() in fortran */
 svdcmp_(q1.f, &i1, &i2, &i3, &i4, q2.f, q3.f);
 
 return 1;
 }
 /*------------------------------------------------------------------------- */
 /* make complex data from a,b: #c=cmake(a,b), 
 	a,b,fltarr & same size. c=a+ib
 */
int ana_cmake(narg, ps)
 int     narg, ps[];
 {
	int ax, ay, i, i0, i1, i2, total_size_in_byte, result_sym;
	int data_type, data_type2, top_type, real_sym, imag_sym;
	int scratch_sym, xq, n1, n2, n;	
	int nd;
	int pd[8];			/* size of dim */
	float *result_ptr, f1, f2;

	real_sym = ps[0];
	/* second arg is optional, if not there, we use symbol 0 (value 0) */
	if(narg >=2) imag_sym = ps[1]; else imag_sym = 0;
	data_type = sym[ps[0]].type;
	if (narg >=2) data_type2 = sym[ps[1]].type; else data_type2 = data_type;
	top_type = MAX(data_type2, data_type);
	
	if (data_type != top_type) { /* someone needs to be converted first */
	 if (data_type < data_type2) xq = real_sym; else xq = imag_sym;
	 switch(top_type)	{
	   case 1: xq = ana_word(1, &xq); break;
	   case 2: xq = ana_long(1, &xq); break;
	   case 3: xq = ana_float(1, &xq); break;
	   case 4: xq = ana_double(1, &xq); break;
	 }
	 if (data_type < data_type2) real_sym = xq; else imag_sym = xq;
	}
	/* everybody is now top_type */
	
	n1 = get_p_c(&real_sym, &q1);
	n2 = get_p_c(&imag_sym, &q2);
	if (n1 == 0 || n2 == 0) return -1;
	/* insist that n1 = n2 */
	n = MAX(n1, n2);
	/* printf("n = %d\n",  n); */
	if (n1 != n2) {
	    printf("CMAKE - mismatch in real/imag input lengths,  %d,  %d\n", 
		n1, n2);
	    return -1;
	}
	/* create the result array, use real_sym unless imag_sym
	 * has more elements
	 */
	if (n1 >= n2) xq = real_sym; else xq = imag_sym;
	switch (sym[xq].class)	{
	case 1:
	case 8:
	    pd[0] = 2;	pd[1]=1; nd = 2; break;
	case 4:
	    h = (struct ahead *) sym[xq].spec.array.ptr;
	    nd = h->ndim;
	    if (nd >= 8) {
		printf("CMAKE - input dimensions must be < 8,  not %d!\n", nd);
		return -1;
	    }
	    for (i=0;i<nd;i++) pd[i+1] = h->dims[i]; pd[0] = 2; nd++;
	    break;
	}

	result_sym = array_scratch(top_type, nd, pd);
	h = (struct ahead *) sym[result_sym].spec.array.ptr;
	q3.l = (int  *) ((char *) h + sizeof( struct ahead));
	switch(top_type)	{
	 case 0: while (n--) {*q3.b++ = *q1.b++; *q3.b++ = *q2.b++; }
		    break;
	 case 1: while (n--) {*q3.w++ = *q1.w++; *q3.w++ = *q2.w++; }
		    break;
	 case 2: while (n--) {*q3.l++ = *q1.l++; *q3.l++ = *q2.l++; }
		    break;
	 case 3: while (n--) {*q3.f++ = *q1.f++; *q3.f++ = *q2.f++; }
		    break;
	 case 4: while (n--) {*q3.d++ = *q1.d++; *q3.d++ = *q2.d++; }
		    break;
	}
	/* printf("end of cmake\n"); */
	return result_sym;	
 }
 /*-------------------------------------------------------------------------*/
int ana_cdiv(narg, ps)
 int     narg, ps[];
 {
        int xsize, ysize=1, result_sym, i;
        int pd[8];                      /* size of dim */
        int ndim, ndim0;                        /* No. of dim of arrays */
        int xstep=1, ystep=1;
        int data_type;
        register float *result_ptr, *cf1, *cf2, *cf3, *cf4;

        data_type = sym[ps[0]].type;
        if(data_type != 3) {
                printf("data_type....%d\n",data_type);
                return 0;               /* not float type  */
        }
        h = (struct ahead *) sym[ps[0]].spec.array.ptr;
        cf1 = (float *) ((char *) h + sizeof( struct ahead));
        cf2 = cf1 + 1;
        ndim = h->ndim;
        if(h->dims[0] != 2) {           /* real data */
                ndim ++;
                xstep = 1;
                xsize = h->dims[0];
 		if (ndim >= 2) ysize = h->dims[1]; else ysize = 1;
       } else {                        /* complex data */
                xstep = 2;
                xsize = h->dims[1];
		if (ndim >= 3) ysize = h->dims[2]; else ysize = 1;
         }

        h = (struct ahead *) sym[ps[1]].spec.array.ptr;
        ndim0 = h->ndim;
        cf3 = (float *) ((char *) h + sizeof( struct ahead));
        cf4 = cf3 + 1;
        if(h->dims[0] != 2) {           /* real data */
                ndim0++;
                ystep = 1;
        } else {
                ystep = 2;
                if(h->ndim > ndim) {
                        ndim = h->ndim;
                        xsize = h->dims[1];
		        if (h->ndim >= 3) ysize = h->dims[2]; else ysize = 1;
                }
        }

        if(ndim < ndim0) ndim = ndim0;
        pd[0] = 2;      pd[1] = xsize;  pd[2] = ysize;
        result_sym = array_scratch(3, ndim, pd);        /* 3: float type */
        h = (struct ahead *) sym[result_sym].spec.array.ptr;
        result_ptr = (float *) ((char *) h + sizeof( struct ahead));

        switch(xstep*10+ystep) {
                case 11:                /* real_x*real_y */
			for(i=0; i<ysize*xsize; i++)    {
                                *result_ptr++ = *cf1++ * *cf3++;
                                *result_ptr++ = 0;
                        }
                        break;
                case 12:                /* real_x*complex_y */
                        for(i=0; i<ysize*xsize; i++) {
				register float fbuf=(*cf3 * *cf3)+(*cf4 * *cf4);
                                *result_ptr++ = (*cf1 * *cf3)/fbuf;
				*result_ptr++ = -(*cf1++ * *cf4)/fbuf;
				cf3 += 2;       cf4 += 2;
			}
			break;
                case 21:                /* complex_x*real_y */
                        for(i=0; i<ysize*xsize; i++) {
				*result_ptr++ = *cf1 / *cf3;
				*result_ptr++ = *cf2 / *cf3++;
				cf1 += 2;       cf2 += 2;
                        }
                        break;
                case 22:                /* complex_x*complex_y */
                        for(i=0; i<ysize*xsize; i++) {
				register float fbuf=(*cf3 * *cf3)+(*cf4 * *cf4);
				*result_ptr++ = 
					((*cf1 * *cf3) + (*cf2 * *cf4) )/fbuf;
				*result_ptr++ = 
					((*cf2 * *cf3) - (*cf1 * *cf4))/fbuf;
				cf1 += 2;       cf2 += 2;
                		cf3 += 2;       cf4 += 2;
        		}
	}
	return result_sym;
	
 }
 /*-------------------------------------------------------------------------*/
int ana_cmul(narg, ps)
 int     narg, ps[];
 {
        int xsize, ysize=1, result_sym, i;
        int pd[8];                      /* size of dim */
        int ndim, ndim0;         		/* No. of dim of arrays */
	int xstep=1, ystep=1;
        int data_type;
        register float *result_ptr, *cf1, *cf2, *cf3, *cf4;

        data_type = sym[ps[0]].type;
	if(data_type != 3) {
		printf("data_type....%d\n",data_type);
		return 0;		/* not float type */
	}
        h = (struct ahead *) sym[ps[0]].spec.array.ptr;
	cf1 = (float *) ((char *) h + sizeof( struct ahead));
        cf2 = cf1 + 1;
	ndim = h->ndim;
	if(h->dims[0] != 2) {		/* real data */
		ndim ++;
		xstep = 1;
		xsize = h->dims[0];
		if (h->ndim >= 2) ysize = h->dims[1]; else ysize = 1;
	} else {			/* complex data */
		xstep = 2;
                xsize = h->dims[1];
		if (h->ndim >= 3) ysize = h->dims[2]; else ysize = 1;
	}

        h = (struct ahead *) sym[ps[1]].spec.array.ptr;
	ndim0 = h->ndim;
	cf3 = (float *) ((char *) h + sizeof( struct ahead));
        cf4 = cf3 + 1;
	if(h->dims[0] != 2) {           /* real data */
		ndim0++;
                ystep = 1;
	} else {    
		ystep = 2;
		if(h->ndim > ndim) {
			ndim = h->ndim;
			xsize = h->dims[1];
			if (h->ndim >= 3) ysize = h->dims[2]; else ysize = 1;
		}
	}

	if(ndim < ndim0) ndim = ndim0;
        pd[0] = 2;      pd[1] = xsize;  pd[2] = ysize;
        result_sym = array_scratch(3, ndim, pd);        /* 3: float type */
	h = (struct ahead *) sym[result_sym].spec.array.ptr;
        result_ptr = (float *) ((char *) h + sizeof( struct ahead));

	switch(xstep*10+ystep) {
		case 11:		/* real_x*real_y */
			for(i=0; i<ysize*xsize; i++)  	{
				*result_ptr++ = *cf1++ * *cf3++;
				*result_ptr++ = 0;
			}
			break;
		case 12:		/* real_x*complex_y */
			for(i=0; i<ysize*xsize; i++) {
				*result_ptr++ = *cf1 * *cf3;
				*result_ptr++ = *cf1++ * *cf4;
				cf3 += 2;       cf4 += 2;
			}
			break;
		case 21:		/* complex_x*real_y */
			for(i=0; i<ysize*xsize; i++) {
				*result_ptr++ = *cf3 * *cf1;
                        	*result_ptr++ = *cf3++ * *cf2;
                        	cf1 += 2;       cf2 += 2;
			}
                        break;
                case 22:		/* complex_x*complex_y */
			for(i=0; i<ysize*xsize; i++) {
				*result_ptr++ = (*cf1 * *cf3) - (*cf2 * *cf4);
				*result_ptr++ = (*cf2 * *cf3) + (*cf1 * *cf4);
				cf1 += 2;       cf2 += 2;
				cf3 += 2;       cf4 += 2;
			}
			break;
	}
			 
        return result_sym;
 }
 /*-------------------------------------------------------------------------*/
int ana_cadd(narg, ps)
 int     narg, ps[];
 {
        int xsize, ysize=1, result_sym, i;
        int pd[8];                      /* size of dim */
        int ndim, ndim0;                        /* No. of dim of arrays */
        int xstep=1, ystep=1;
        int data_type;
        register float *result_ptr, *cf1, *cf2, *cf3, *cf4;

        data_type = sym[ps[0]].type;
        if(data_type != 3) {
                printf("data_type....%d\n",data_type);
                return 0;               /* not float type */
        }
        h = (struct ahead *) sym[ps[0]].spec.array.ptr;
        cf1 = (float *) ((char *) h + sizeof( struct ahead));
        cf2 = cf1 + 1;
        ndim = h->ndim;
        if(h->dims[0] != 2) {           /* real data */
                ndim ++;
                xstep = 1;
                xsize = h->dims[0];
		if (h->ndim >= 2) ysize = h->dims[1]; else ysize = 1;
        } else {                        /* complex data */
                xstep = 2;
                xsize = h->dims[1];
		if (ndim >= 3) ysize = h->dims[2]; else ysize = 1;
        }

        h = (struct ahead *) sym[ps[1]].spec.array.ptr;
        ndim0 = h->ndim;
        cf3 = (float *) ((char *) h + sizeof( struct ahead));
        cf4 = cf3 + 1;
        if(h->dims[0] != 2) {           /* real data */
                ndim0++;
                ystep = 1;
        } else {
                ystep = 2;
                if(h->ndim > ndim) {
                        ndim = h->ndim;
                        xsize = h->dims[1];
			if (h->ndim >= 3) ysize = h->dims[2]; else ysize = 1;

                }
        }

        if(ndim < ndim0) ndim = ndim0;
        pd[0] = 2;      pd[1] = xsize;  pd[2] = ysize;
        result_sym = array_scratch(3, ndim, pd);        /* 3: float type */
        h = (struct ahead *) sym[result_sym].spec.array.ptr;
        result_ptr = (float *) ((char *) h + sizeof( struct ahead));

        /* printf("xstep*10+ystep = %d\n",  xstep*10+ystep); */
	switch(xstep*10+ystep) {
                case 11:                /* real_x*real_y */
                        for(i=0; i<ysize*xsize; i++)    {
				*result_ptr++ = *cf1++ + *cf3++;
                                *result_ptr++ = 0;
                        }
                        break;
                case 12:                /* real_x*complex_y */
                        for(i=0; i<ysize*xsize; i++) {
                                *result_ptr++ = *cf1++ +  *cf3;
                                *result_ptr++ = *cf4;
                                cf3 += 2;       cf4 += 2;
                        }
                        break;
                case 21:                /* complex_x*real_y */
                        for(i=0; i<ysize*xsize; i++) {
                                *result_ptr++ = *cf1 + *cf3++;
				*result_ptr++ = *cf2;
                                cf1 += 2;       cf2 += 2;
                        }
                        break;
                case 22:                /* complex_x*complex_y */
                        for(i=0; i<ysize*xsize; i++) {
				*result_ptr++ = *cf1 + *cf3;
				*result_ptr++ = *cf2 + *cf4;
				cf1 += 2;       cf2 += 2;
                                cf3 += 2;       cf4 += 2;
                        }
                        break;
        }

        return result_sym;
 }
 /*-------------------------------------------------------------------------*/
int ana_csub(narg, ps)		/* complex sub */
 int     narg, ps[];
 {
        int xsize, ysize=1, result_sym, i;
        int pd[8];                      /* size of dim */
        int ndim, ndim0;                        /* No. of dim of arrays */
        int xstep=1, ystep=1;
        int data_type;
        register float *result_ptr, *cf1, *cf2, *cf3, *cf4;

        data_type = sym[ps[0]].type;
        if(data_type != 3) {
                printf("data_type....%d\n",data_type);
                return 0;               /* not float type */
        }
        h = (struct ahead *) sym[ps[0]].spec.array.ptr;
        cf1 = (float *) ((char *) h + sizeof( struct ahead));
        cf2 = cf1 + 1;
        ndim = h->ndim;
        if(h->dims[0] != 2) {           /* real data */
                ndim ++;
                xstep = 1;
                xsize = h->dims[0];
		if (h->ndim >= 2) ysize = h->dims[1]; else ysize = 1;
        } else {                        /* complex data */
                xstep = 2;
                xsize = h->dims[1];
        	if (h->ndim >= 3) ysize = h->dims[2]; else ysize = 1;
    }

        h = (struct ahead *) sym[ps[1]].spec.array.ptr;
        ndim0 = h->ndim;
        cf3 = (float *) ((char *) h + sizeof( struct ahead));
        cf4 = cf3 + 1;
        if(h->dims[0] != 2) {           /* real data */
                ndim0++;
                ystep = 1;
        } else {
                ystep = 2;
                if(h->ndim > ndim) {
                        ndim = h->ndim;
                        xsize = h->dims[1];
        		if (h->ndim >= 3) ysize = h->dims[2]; else ysize = 1;
                }
        }

                                                        
        pd[0] = 2;      pd[1] = xsize;  pd[2] = ysize;
        result_sym = array_scratch(3, ndim, pd);        /* 3: float type */
        h = (struct ahead *) sym[result_sym].spec.array.ptr;
        result_ptr = (float *) ((char *) h + sizeof( struct ahead));

        switch(xstep*10+ystep) {
                case 11:                /* real_x*real_y */
                        for(i=0; i<ysize*xsize; i++)    {
                                *result_ptr++ = *cf1++ - *cf3++;
                                *result_ptr++ = 0;
                        }
                        break;
                case 12:                /* real_x*complex_y */
                        for(i=0; i<ysize*xsize; i++) {
                                *result_ptr++ = *cf1++ - *cf3;
                                *result_ptr++ = *cf4;
                                cf3 += 2;       cf4 += 2;
                        }
                        break;
                case 21:                /* complex_x*real_y */
                        for(i=0; i<ysize*xsize; i++) {
                                *result_ptr++ = *cf1 - *cf3++;
                                *result_ptr++ = *cf2;
                                cf1 += 2;       cf2 += 2;
                        }
                        break;
                case 22:                /* complex_x*complex_y */
                        for(i=0; i<ysize*xsize; i++) {
                                *result_ptr++ = *cf1 - *cf3;
                                *result_ptr++ = *cf2 - *cf4;
                                cf1 += 2;       cf2 += 2;
                                cf3 += 2;       cf4 += 2;
                        }
                        break;
        }

        return result_sym;
}
 /*-------------------------------------------------------------------------*/
int ana_cconj(narg, ps)		/* complex conjugate a+jb -> a-jb */
 int     narg, ps[];
 {
        int i, ndim, xsize, ysize;
	register float *cf1;

	h = (struct ahead *) sym[ps[0]].spec.array.ptr;
	if(h->dims[0] != 2) {		/* real data, just return */
		printf("\007Cann't do conjugate with real data\n");
		return ps[0];
	} else {
		cf1 = (float *) ((char *) h + sizeof( struct ahead)) + 1;
		ndim = h->ndim;
		xsize = h->dims[1];
		if (ndim >= 3) ysize = h->dims[2]; else ysize = 1;
		if(ysize == 0 || ndim == 1) ysize = 1;
		for(i=0; i<ysize*xsize; i++) 	{
			*cf1 = - *cf1;
			cf1 += 2;
		}
	}
	return ps[0];
 }
 /*-------------------------------------------------------------------------*/
int ana_csqr(narg, ps)
 int     narg, ps[];
 {
        int xsize, ysize=1, result_sym, i;
        int pd[8];                      /* size of dim */
        int ndim;                       /* No. of dim of arrays */
        int xstep=1, ystep=1;
        register float *result_ptr, *cf1, *cf2;

        h = (struct ahead *) sym[ps[0]].spec.array.ptr;
        cf1 = (float *) ((char *) h + sizeof( struct ahead));
        cf2 = cf1 + 1;
        ndim = h->ndim;
        if(h->dims[0] != 2) {           /* real data */
                xstep = 1;
                xsize = h->dims[0]; ysize = h->dims[1];
                if(ysize == 0) ysize = 1;
        } else {                        /* complex data */
                xstep = 2;
                xsize = h->dims[1];
		if (ndim >= 3) ysize = h->dims[2]; else ysize = 1;
                if(ysize == 0 ) ysize = 1;
		ndim--;
        }

	pd[0] = xsize;  pd[1] = ysize;
        result_sym = array_scratch(3, ndim, pd);        /* 3: float type */
        h = (struct ahead *) sym[result_sym].spec.array.ptr;
        result_ptr = (float *) ((char *) h + sizeof( struct ahead));

	switch(xstep) {
                case 1:                /* real_x */
                        for(i=0; i<ysize*xsize; i++) {
                                *result_ptr++ = *cf1 * *cf1;
				cf1++;
			}
			break;
		case 2:			/*complex_x */
			for(i=0; i<ysize*xsize; i++)	{
				*result_ptr++ = (*cf1 * *cf1) + (*cf2 * *cf2);
				cf1 += 2;       cf2 += 2;
			}
			break;
	}

	return result_sym;
 }