mathface.cpp

普林斯顿开发的快速球面调和变换算法

// MathFace.cpp: implementation of the MathFace class.
//
//////////////////////////////////////////////////////////////////////


#include "stdafx.h"
#include "MathFace.h"

#include "stdio.h"

extern "C"
{
	int getc( FILE *stream );
	int isdigit( int c );
	int ungetc( int c, FILE *stream );
    int fscanf( FILE *stream, const char *format,... );


}

#define LOG_10_2 0.3010299956639812  /* log base 10 of 2 */

int MMChopLevel = -15;  /* numbers smaller than 10^ChopLevel are treated as 0 */
int MMPrecision = 16;    /* number of digits after the radix */



#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
extern "C"
{

	int putc( int c, FILE *stream );
	
}


MathFace::MathFace()
{

}

MathFace::~MathFace()
{

}

int MathFace::readNumber(FILE *f, double *real, double *imag)
{
	  int c, gotImag, foo, scale;
  double negate;

  fscanf(f, " ");

  gotImag = 0;
  if((c = getc(f)) == '-') {
    negate = -1.0;
    fscanf(f, " ");
  } else {
    negate = 1.0;
    ungetc(c, f);
  }


  /* added by Sean */
  if ((c = getc(f)) == '(')
    fscanf(f," ");
  else
    ungetc(c,f);
  /* end of added by Sean */

  if((c = getc(f)) == 'I') {
    *real = 0.0;
    *imag = negate;
    gotImag = 1;
    fscanf(f, " ");
  } else if((c == '.') || isdigit(c)) {
    ungetc(c, f);
    foo = fscanf(f, " %lf ", real);
    *real *= negate;
    /*
     * scanf only reads scientific notation with small e's, so I have
     *  to catch capital E's
     */
    if(((c = getc(f)) != '*') && (c != 'E')) {
      ungetc(c, f);
    } else {
      fscanf(f, " ");
      if(c == 'E') {
	if((c = getc(f)) != '+') {
	  ungetc(c, f);
	}
	fscanf(f, "%d ", &scale);
	*real *= pow(10.0, (double) scale);	
      } else if((c = getc(f)) == '1') {
	fscanf(f, "0 ^ ");
	if((c = getc(f)) != '+') {
	  ungetc(c, f);
	}
	fscanf(f, "%d ", &scale);
	*real *= pow(10.0, (double) scale);

	/* added by Sean */
	if ((c = getc(f)) == ')')
	  fscanf(f," ");
	else
	  ungetc(c,f);
	/* end of added by Sean */

      } else if(c == 'I') {
	*imag = *real;
	*real = 0.0;
	gotImag = 1;
	fscanf(f, " ");
      } else {
	ungetc(c, f);
      }
    }
  } else {
    ungetc(c, f);
    return -1;
  }

  if((((c = getc(f)) == '+') || (c == '-')) && (gotImag == 0)) {
    if(c == '-') {
      negate = -1.0;
    } else {
      negate = 1.0;
    }
    fscanf(f, " ");
    if((c = getc(f)) == 'I') {
      *imag = negate;
      fscanf(f, " ");
    } else if((c == '.') || isdigit(c)) {
      ungetc(c, f);
      fscanf(f, " %lf ", imag);
      *imag *= negate;
      /*
       * scanf only reads scientific notation with small e's, so I have
       *  to catch capital E's
       */
      if((c = getc(f)) == 'E') {
	fscanf(f, " ");
	if((c = getc(f)) != '+') {
	  ungetc(c, f);
	}
	fscanf(f, "%d * ", &scale);
	*imag *= pow(10.0, (double) scale);
      } else if(c != '*') {
	return -1;
      }
      fscanf(f, " ");
      if((c = getc(f)) != '1') {
	ungetc(c, f);
      } else {
	fscanf(f, "0 ^ ");
	if((c = getc(f)) != '+') {
	  ungetc(c, f);
	}
	fscanf(f, "%d * ", &scale);
	*imag *= pow(10.0, (double) scale);
      }
      if((c = getc(f)) != 'I') {
	ungetc(c, f);
      }
      fscanf(f, " ");
    } else {
      ungetc(c, f);
      return -1;
    }
  } else if((c == '*') && (gotImag == 0)) {
    fscanf(f, " ");
    if((c = getc(f)) == 'I') {
      *imag = *real;
      *real = 0.0;
      fscanf(f, " ");
    } else {
      ungetc(c, f);
    }
  } else {
    ungetc(c, f);
    if(gotImag == 0) {
      *imag = 0.0;
    }
  }

  return 0;

}

void MathFace::printNumber(FILE *f, double real, double imag)
{


}

void MathFace::seanprintNumber(FILE *f, double real, double imag)
{
	fprintf(f,"%17.15f + I * %17.15f",real,imag);

}

 /*
 读文件,如果表是1-D,那么高度为0
 返回值:1,没有错误;0,有问题
 */
int MathFace::readMMRealTable(FILE *f, double *table, int size, int *width, int *height)
{
	  int i,j, foo, n, c;
	  int flat;
  double real, imag;

  n = 0;
  foo = fscanf(f, " { ");
  //foo=getc(f);
//fscanf(f," ",&foo);
  
  //准备接收第一个字符
  if((c = getc(f)) != '{') {
    ungetc(c, f);
    flat = 1;
  } else {
    fscanf(f, " ");
    flat = 0;
  }
  j = 0;

  while(((c = getc(f)) != '}') && (j < size)) {
	ungetc(c, f);
    readNumber(f, &real, &imag);
    if(imag != 0.0) {
      return 0;
    }
    table[n] = real;
    n++;
    j++;
    if((c = getc(f)) != ',') {
      ungetc(c, f);
    } else {
      fscanf(f, " ");
    }
  }
  if((j == size) && (c != '}')) { return 0; }
  *width = j;
  if(flat == 0) {
    fscanf(f, " ");
    i = 1;
    while(((c = getc(f)) != '}') && (i*(*width) < size)) {
      foo = fscanf(f, " { ");
      j = 0;
      while(((c = getc(f)) != '}') && (j < *width)) {
	ungetc(c, f);
	readNumber(f, &real, &imag);
	if(imag != 0.0) {
	  return 0;
	}
	table[n] = real;
	n++;
	j++;
	if((c = getc(f)) != ',') {
	  ungetc(c, f);
	} else {
	  fscanf(f, " ");
	}
      }
      if(c != '}') {
	return 0;
      }
      fscanf(f, " ");
      i++;
    }
    if((i*(*width) == size) && (c != '}')) { return 0; }
    *height = i;
  } else {
    *height = 0;
  }

  return 1;


}

void MathFace::printMMRealTable(FILE *f, double *table, int width, int height)
{

}

int MathFace::readMMComplexTable(FILE *f, double *real, double *imag, int size, int *width, int *height, int *isReal)
{
	return 0;

}


void MathFace::printMMComplexTable(FILE *f, double *real, double *imag, int width, int height)
{
	

}

void MathFace::seanprintMMComplexTable(FILE *f, double *real, double *imag, int width, int height)
{

}

void MathFace::seanprintMMRealTable(FILE *f, double *table, int width, int height)
{
	  int i, j, flat;


  if(height != 0) {
    fprintf(f, "{");
    flat = 0;
  } else {
    height = 1;
    flat = 1;
  }

  for(i = 0; i < height; i++) {
    if(i != 0) {
      fprintf(f, ",\n ");
    }
    fprintf(f, "{");
    for(j = 0; j < width; j++) {
      if(j != 0) {
	putc(',', f);
      }
      seanprintNumber(f, table[width*i+j], 0.0);
    }
    fprintf(f, "}");
    fflush(f);
  }
  if(flat == 0) {
    fprintf(f, "}\n");
  } else {
    putc('\n', f);
  }


}