class CAtom{
public:
	int    number;
	double charge;
	ogVector position;
	void loadCube( FILE *fp ){ fscanf( fp, "%d%lf%lf%lf%lf", &number, &charge, &position.x, &position.y, &position.z ); }
	void saveCube( FILE *fp ){ fprintf( fp, "%d %f %f %f %f\n", number, float(charge), float(position.x), float(position.y), float(position.z) ); }
	void loadBin( FILE *fp ){
		float c = charge; float x = position.x; float y = position.y; float z = position.z;
		fread( (int   *)&number, 1, sizeof(int), fp );
		fread( (float *)&c, 1, sizeof(float), fp );
		fread( (float *)&x, 1, sizeof(float), fp );
		fread( (float *)&y, 1, sizeof(float), fp );
		fread( (float *)&z, 1, sizeof(float), fp );
	}
	void saveBin( FILE *fp ){
		float c = charge; float x = position.x; float y = position.y; float z = position.z;
		fwrite( (int   *)&number, 1, sizeof(int), fp );
		fwrite( (float *)&c, 1, sizeof(float), fp );
		fwrite( (float *)&x, 1, sizeof(float), fp );
		fwrite( (float *)&y, 1, sizeof(float), fp );
		fwrite( (float *)&z, 1, sizeof(float), fp );
	}
};
//-------------------------------------------------------------------------------------------
class CCube{
public:
	int     IFlag;
	ogVector origin;
	int     xSize;
	ogVector xUnit;
	int     ySize;
	ogVector yUnit;
	int     zSize;
	ogVector zUnit;
	ogVector center;
	CAtom  *atom;
	int     num1;
	int     num2;
	float  *d;
	GLuint *texName;

	float density( int x, int y, int z ){ return d[z + (y + x*ySize) * zSize]; }
	void convert( double k1, double k2 ){
		for( long i=0; i<xSize*ySize*zSize; i++ ) d[i] = (d[i] > 0 ? 1.0 : -1.0) * k1 * exp( k2 * log( fabsf(d[i]) ) );
	}
	void construct(){ atom = NULL; d = NULL; texName = NULL; }
	void construct( char *filename, double k1, double k2 ){ this->construct(); this->load( filename ); this->convert( k1, k2 ); }
	CCube(){ this->construct(); }
	CCube( char *filename, double k1, double k2 ){ this->construct( filename, k1, k2 ); }
	~CCube(){ if( atom != NULL ) delete atom; if( d != NULL ) delete d; if( texName != NULL ) delete texName; }
	void glTexImage( ogImage &image, int textureNumber );
	void texImage();

	void drawYZ( int x );
	void drawZX( int y );
	void drawXY( int z );
	void drawYZ(){ for( int x=0; x<xSize; x++ ){ this->drawYZ( x ); }}
	void drawZY(){ for( int x=xSize-1; x>=0; x-- ){ this->drawYZ( x ); }}
	void drawZX(){ for( int y=0; y<ySize; y++ ){ this->drawZX( y ); }}
	void drawXZ(){ for( int y=ySize-1; y>=0; y-- ){ this->drawZX( y ); }}
	void drawXY(){ for( int z=0; z<zSize; z++ ){ this->drawXY( z ); }}
	void drawYX(){ for( int z=zSize-1; z>=0; z-- ){ this->drawXY( z ); }}

	void drawYZ( int x, double intensity, double ratio, CCube &cube );
	void drawZX( int y, double intensity, double ratio, CCube &cube );
	void drawXY( int z, double intensity, double ratio, CCube &cube );
	void drawYZ( double intensity, double ratio, CCube &cube ){ for( int x=0; x<xSize; x++ ) this->drawYZ( x, intensity, ratio, cube ); }
	void drawZY( double intensity, double ratio, CCube &cube ){ for( int x=xSize-1; x>=0; x-- ) this->drawYZ( x, intensity, ratio, cube ); }
	void drawZX( double intensity, double ratio, CCube &cube ){ for( int y=0; y<ySize; y++ ) this->drawZX( y, intensity, ratio, cube ); }
	void drawXZ( double intensity, double ratio, CCube &cube ){ for( int y=ySize-1; y>=0; y-- ) this->drawZX( y, intensity, ratio, cube ); }
	void drawXY( double intensity, double ratio, CCube &cube ){ for( int z=0; z<zSize; z++ ) this->drawXY( z, intensity, ratio, cube ); }
	void drawYX( double intensity, double ratio, CCube &cube ){ for( int z=zSize-1; z>=0; z-- ) this->drawXY( z, intensity, ratio, cube ); }
	void draw( ogVector &direction, double intensity );
	void draw( ogVector &direction, double intensity, double ratio, CCube &cube );

	void drawAtoms();
	void drawAtoms( double ratio, CCube &cube );
	void drawBonds();
	void drawBonds( double ratio, CCube &cube );
	void drawFrame( GLubyte red, GLubyte green, GLubyte blue );

	int loadCube( char *filename );
	int saveCube( char *filename );
	int loadBin( char *filename );
	int saveBin( char *filename );
	int load( char *filename );
	int save( char *filename );

	void interpolate( CCube &cube1, CCube &cube2, double a0, double b0 ){
		double a = a0 / (a0+b0);
		double b = b0 / (a0+b0);
		for( long i=0; i<xSize*ySize*zSize; i++ ) d[i] = cube1.d[i] * a + cube2.d[i] * b;
	}
};

void CCube::glTexImage( ogImage &image, int textureNumber ){
	glGenTextures( 1, &texName[textureNumber] );
	glBindTexture( GL_TEXTURE_2D, texName[textureNumber] );
	glPixelStorei( GL_UNPACK_ALIGNMENT,1 );
	glTexParameterf( GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT );
	glTexParameterf( GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT );
	glTexParameterf( GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR );
	glTexParameterf( GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR );
	glTexEnvf( GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE );
	glTexImage2D( GL_TEXTURE_2D, 0, image.format(), image.width, image.height, 0, image.format(), GL_UNSIGNED_BYTE, image.array );
}

void CCube::texImage(){
	if( texName != NULL ) return;
	texName = new GLuint[xSize+ySize+zSize];

	ogImage image( ySize, zSize, 4 );
	ogImage image128( 128, 128, 4 );
	for( int x=0; x<xSize; x++ ){
		for( int z=0; z<zSize; z++ ){
		for( int y=0; y<ySize; y++ ){
			if( this->density( x, y, z ) > 0 ) image.set( y, z, 255, 255, 255, cramp( int(this->density( x, y, z )) ) );
			else image.set( y, z, 255, 0, 0, cramp( int(-this->density( x, y, z )) ) );
		}}
		gluScaleImage( GL_RGBA, ySize, zSize, GL_UNSIGNED_BYTE, image.array, 128, 128, GL_UNSIGNED_BYTE, image128.array );
		this->glTexImage( image128, x );
	}
	image.construct( zSize, xSize, 4 );
	for( int y=0; y<ySize; y++ ){
		for( int x=0; x<xSize; x++ ){
		for( int z=0; z<zSize; z++ ){
			if( this->density( x, y, z ) > 0 ) image.set( z, x, 255, 255, 255, cramp( int(this->density( x, y, z )) ) );
			else image.set( z, x, 255, 0, 0, cramp( int(-this->density( x, y, z )) ) );
		}}
		gluScaleImage( GL_RGBA, zSize, xSize, GL_UNSIGNED_BYTE, image.array, 128, 128, GL_UNSIGNED_BYTE, image128.array );
		this->glTexImage( image128, xSize + y );
	}
	image.construct( xSize, ySize, 4 );
	for( int z=0; z<zSize; z++ ){
		for( int y=0; y<ySize; y++ ){
		for( int x=0; x<xSize; x++ ){
			if( this->density( x, y, z ) > 0 ) image.set( x, y, 255, 255, 255, cramp( int(this->density( x, y, z )) ) );
			else image.set( x, y, 255, 0, 0, cramp( int(-this->density( x, y, z )) ) );
		}}
		gluScaleImage( GL_RGBA, xSize, ySize, GL_UNSIGNED_BYTE, image.array, 128, 128, GL_UNSIGNED_BYTE, image128.array );
		this->glTexImage( image128, xSize + ySize + z );
	}

	if( d != NULL ) delete d;
	d = NULL;
}

void CCube::drawYZ( int x ){
	glBindTexture( GL_TEXTURE_2D, texName[x] );
	glBegin(GL_POLYGON);
	glTexCoord2i( 0, 0 );
	glVertex3d( double(x) * xUnit.abs(), 0, 0 );
	glTexCoord2i( 0, 1 );
	glVertex3d( double(x) * xUnit.abs(), 0, zSize * zUnit.abs() );
	glTexCoord2i( 1, 1 );
	glVertex3d( double(x) * xUnit.abs(), ySize * yUnit.abs(), zSize * zUnit.abs() );
	glTexCoord2i( 1, 0 );
	glVertex3d( double(x) * xUnit.abs(), ySize * yUnit.abs(), 0 );
	glEnd();
}

void CCube::drawZX( int y ){
	glBindTexture( GL_TEXTURE_2D, texName[xSize+y] );
	glBegin(GL_POLYGON);
	glTexCoord2i( 0, 0 );
	glVertex3d( 0, double( y ) * yUnit.abs(), 0 );
	glTexCoord2i( 0, 1 );
	glVertex3d( xSize * xUnit.abs(), double( y ) * yUnit.abs(), 0 );
	glTexCoord2i( 1, 1 );
	glVertex3d( xSize * xUnit.abs(), double( y ) * yUnit.abs(), zSize * zUnit.abs());
	glTexCoord2i( 1, 0 );
	glVertex3d( 0, double( y ) * yUnit.abs(), zSize * zUnit.abs());
	glEnd();
}

void CCube::drawXY( int z ){
	glBindTexture( GL_TEXTURE_2D, texName[xSize+ySize+z] );
	glBegin(GL_POLYGON);
	glTexCoord2i( 0, 0 );
	glVertex3d( 0, 0, double( z ) * zUnit.abs());
	glTexCoord2i( 0, 1 );
	glVertex3d( 0, ySize * yUnit.abs(), double( z ) * zUnit.abs() );
	glTexCoord2i( 1, 1 );
	glVertex3d( xSize * xUnit.abs(), ySize * yUnit.abs(), double( z ) * zUnit.abs() );
	glTexCoord2i( 1, 0 );
	glVertex3d( xSize * xUnit.abs(), 0, double( z ) * zUnit.abs() );
	glEnd();
}

void CCube::drawYZ( int x, double intensity, double ratio, CCube &cube ){
	glBlendFunc( GL_ZERO,GL_ONE_MINUS_SRC_ALPHA );
	glColor4ub( 255, 255, 255, cramp(int(255.0*intensity)) );
	ratio < 0.5 ? this->drawYZ( x ) : cube.drawYZ( x );
	glBlendFunc( GL_SRC_ALPHA,GL_ONE );
	glColor4ub( 255, 255, 255, cramp(int(255.0*intensity*(1.0-ratio))) ); this->drawYZ( x );
	glBlendFunc( GL_SRC_ALPHA,GL_ONE );
	glColor4ub( 255, 255, 255, cramp(int(255.0*intensity*ratio)) ); cube.drawYZ( x );
}

void CCube::drawZX( int y, double intensity, double ratio, CCube &cube ){
	glBlendFunc( GL_ZERO,GL_ONE_MINUS_SRC_ALPHA );
	glColor4ub( 255, 255, 255, cramp(int(255.0*intensity)) );
	ratio < 0.5 ? this->drawZX( y ) : cube.drawZX( y );
	glBlendFunc( GL_SRC_ALPHA,GL_ONE );
	glColor4ub( 255, 255, 255, cramp(int(255.0*intensity*(1.0-ratio))) ); this->drawZX( y );
	glBlendFunc( GL_SRC_ALPHA,GL_ONE );
	glColor4ub( 255, 255, 255, cramp(int(255.0*intensity*ratio)) ); cube.drawZX( y );
}

void CCube::drawXY( int z, double intensity, double ratio, CCube &cube ){
	glBlendFunc( GL_ZERO,GL_ONE_MINUS_SRC_ALPHA );
	glColor4ub( 255, 255, 255, cramp(int(255.0*intensity)) );
	ratio < 0.5 ? this->drawXY( z ) : cube.drawXY( z );
	glBlendFunc( GL_SRC_ALPHA,GL_ONE );
	glColor4ub( 255, 255, 255, cramp(int(255.0*intensity*(1.0-ratio))) ); this->drawXY( z );
	glBlendFunc( GL_SRC_ALPHA,GL_ONE );
	glColor4ub( 255, 255, 255, cramp(int(255.0*intensity*ratio)) ); cube.drawXY( z );
}

void CCube::draw( ogVector &direction, double intensity ){
	if( texName == NULL ) this->texImage();
	glBlendFunc( GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA );
	glColor4ub( 255, 255, 255, cramp(int(255*intensity)) );
	og.PushMatrix();
	og.Translate(-center);
	if( fabs(direction.x) > fabs(direction.y) && fabs(direction.x) > fabs(direction.z) ){
		if( direction.x < 0 ) this->drawYZ();
		else                  this->drawZY();
	}else if( fabs(direction.y) > fabs(direction.z) ){
		if( direction.y < 0 ) this->drawZX();
		else                  this->drawXZ();
	}else{
		if( direction.z < 0 ) this->drawXY();
		else                  this->drawYX();
	}
	og.PopMatrix();
}

void CCube::draw( ogVector &direction, double intensity, double ratio, CCube &cube ){
	if( texName == NULL ) this->texImage();
	if( cube.texName == NULL ) cube.texImage();
	og.PushMatrix();
	og.Translate(-center);
	if( fabs(direction.x) > fabs(direction.y) && fabs(direction.x) > fabs(direction.z) ){
		if( direction.x < 0 ) this->drawYZ( intensity, ratio, cube );
		else                  this->drawZY( intensity, ratio, cube );
	}else if( fabs(direction.y) > fabs(direction.z) ){
		if( direction.y < 0 ) this->drawZX( intensity, ratio, cube );
		else                  this->drawXZ( intensity, ratio, cube );
	}else{
		if( direction.z < 0 ) this->drawXY( intensity, ratio, cube );
		else                  this->drawYX( intensity, ratio, cube );
	}
	og.PopMatrix();
}

void CCube::drawAtoms(){
	og.PushMatrix();
	og.Translate(-(origin+center));
	for(int i=0; i<abs(IFlag); i++){
		ogVector position = atom[i].position;
		switch( atom[i].number ){
			case  1:og.DrawSphere(position, 0.3, 255,255,255, 2); break;
			case  6:og.DrawSphere(position, 0.3, 255,0,0, 2); break;
			default:og.DrawSphere(position, 0.3, 0,0,255, 2);; break;
		}
	}
	og.PopMatrix();
}

void CCube::drawAtoms( double ratio, CCube &cube ){
	og.PushMatrix();
	og.Translate(-(origin+center));
	for(int i=0; i<abs(IFlag); i++){
		ogVector position = atom[i].position*(1.0-ratio) + cube.atom[i].position*ratio;
		switch( atom[i].number ){
			case  1:og.DrawSphere(position, 0.3, 255,255,255, 2); break;
			case  6:og.DrawSphere(position, 0.3, 255,0,0, 2); break;
			default:og.DrawSphere(position, 0.3, 0,0,255, 2);; break;
		}
	}
	og.PopMatrix();
}

void CCube::drawBonds(){
	double radius[9] = { 0.0, 0.03, 0.0, 0.0, 0.0, 0.0, 0.077, 0.074, 0.074 };
	for( int i=0; i<9; i++ ) radius[i] /= (0.0523167+0.007);

	og.PushMatrix();
	og.Translate(-(origin+center));
	for( i=0; i<abs(IFlag); i++ ){	
		ogVector p1 = atom[i].position;
		for( int j=i+1; j<abs(IFlag); j++){
			ogVector p2 = atom[j].position;
			if( radius[atom[i].number] + radius[atom[j].number] + 0.5 > (p2-p1).abs() ){
				og.DrawCylinder(p1, p2, 0.1, 200, 200, 200, 16);
			}
		}
	}
	og.PopMatrix();
}

void CCube::drawBonds( double ratio, CCube &cube ){
	double radius[9] = { 0.0, 0.03, 0.0, 0.0, 0.0, 0.0, 0.077, 0.074, 0.074 };
	for( int i=0; i<9; i++ ) radius[i] /= (0.0523167+0.007);

	og.PushMatrix();
	og.Translate(-(origin+center));
	for( i=0; i<abs(IFlag); i++ ){	
		for( int j=i+1; j<abs(IFlag); j++){
			ogVector p1 = atom[i].position; // p1 must be refreshed every time
			ogVector p2 = atom[j].position;
			if( radius[atom[i].number] + radius[atom[j].number] + 0.5 > (p2-p1).abs() ){
				p1 *= 1.0-ratio;
				p1 += cube.atom[i].position*ratio;
				p2 *= 1.0-ratio;
				p2 += cube.atom[j].position*ratio;
				//p1 -= origin + center;
				//p2 -= origin + center;
				og.DrawCylinder(p1, p2, 0.1, 200, 200, 200, 16);
			}
		}
	}
	og.PopMatrix();
}

void CCube::drawFrame( GLubyte red, GLubyte green, GLubyte blue ){
	ogVector v[8];
	v[0] = ( - xUnit * double(xSize-1)/2.0 - yUnit * double(ySize-1)/2.0 - zUnit * double(zSize-1)/2.0 );
	v[1] = ( - xUnit * double(xSize-1)/2.0 + yUnit * double(ySize-1)/2.0 - zUnit * double(zSize-1)/2.0 );
	v[2] = (   xUnit * double(xSize-1)/2.0 + yUnit * double(ySize-1)/2.0 - zUnit * double(zSize-1)/2.0 );
	v[3] = (   xUnit * double(xSize-1)/2.0 - yUnit * double(ySize-1)/2.0 - zUnit * double(zSize-1)/2.0 );
	v[4] = ( - xUnit * double(xSize-1)/2.0 - yUnit * double(ySize-1)/2.0 + zUnit * double(zSize-1)/2.0 );
	v[5] = ( - xUnit * double(xSize-1)/2.0 + yUnit * double(ySize-1)/2.0 + zUnit * double(zSize-1)/2.0 );
	v[6] = (   xUnit * double(xSize-1)/2.0 + yUnit * double(ySize-1)/2.0 + zUnit * double(zSize-1)/2.0 );
	v[7] = (   xUnit * double(xSize-1)/2.0 - yUnit * double(ySize-1)/2.0 + zUnit * double(zSize-1)/2.0 );
	for(int i=0; i<4; i++){
		og.DrawLine(v[i],v[(i+1)%4],red,green,blue);
		og.DrawLine(v[i+4],v[(i+1)%4+4],red,green,blue);
		og.DrawLine(v[i],v[i+4],red,green,blue);
	}
}

int CCube::loadCube( char *filename ){
FILE *fp;
char ch;
	if( (fp = fopen( filename, "r" )) == NULL ){
		printf("%s not found!n", filename);
		return 1;
	}
	while( (ch=fgetc(fp)) != '\n' );
	while( (ch=fgetc(fp)) != '\n' );

	fscanf( fp, "%d %lf %lf %lf", &IFlag, &origin.x, &origin.y, &origin.z );
	fscanf( fp, "%d %lf %lf %lf", &xSize, &xUnit.x, &xUnit.y, &xUnit.z );
	fscanf( fp, "%d %lf %lf %lf", &ySize, &yUnit.x, &yUnit.y, &yUnit.z );
	fscanf( fp, "%d %lf %lf %lf", &zSize, &zUnit.x, &zUnit.y, &zUnit.z );
	if( atom != NULL ) delete atom;
	atom = new CAtom[abs(IFlag)];
	for( long i=0; i< abs(IFlag); i++ ) atom[i].loadCube( fp );
	if( IFlag < 0 ) fscanf( fp, "%d %d", &num1, &num2 );
	if( d != NULL ) delete d;
	d = new float[xSize*ySize*zSize];
	for( i=0; i<xSize*ySize*zSize; i++ ) fscanf( fp, "%e", &d[i] );
	fclose( fp );
	printf( "%s loaded\n", filename );
	return 0;
}

int CCube::saveCube( char *filename ){
FILE *fp;
	if( (fp = fopen( filename, "w" )) == NULL ) return 1;
	fprintf( fp, "\n" );
	fprintf( fp, "\n" );

	fprintf( fp, "%d %lf %lf %lf\n", IFlag, origin.x, origin.y, origin.z );
	fprintf( fp, "%d %lf %lf %lf\n", xSize, xUnit.x, xUnit.y, xUnit.z );
	fprintf( fp, "%d %lf %lf %lf\n", ySize, yUnit.x, yUnit.y, yUnit.z );
	fprintf( fp, "%d %lf %lf %lf\n", zSize, zUnit.x, zUnit.y, zUnit.z );
	for( long i=0; i< abs(IFlag); i++ ) atom[i].saveCube( fp );
	if( IFlag < 0 ) fprintf( fp, "%d %d\n", num1, num2 );
	for( i=0; i< xSize*ySize*zSize; i++ ) fprintf( fp, "%e ", d[i] );
	fclose( fp );
	return 0;
}

int CCube::loadBin( char *filename ){
FILE *fp;
	if( (fp = fopen( filename, "rb" )) == NULL ){
		printf("%s not found!n", filename);
		return 1;
	}
	fread( (int *)&IFlag, 1, sizeof(int), fp );
	fread( (double *)&origin.x, 1, sizeof(double), fp );
	fread( (double *)&origin.y, 1, sizeof(double), fp );
	fread( (double *)&origin.z, 1, sizeof(double), fp );
	fread( (int *)&xSize, 1, sizeof(int), fp );
	fread( (double *)&xUnit.x, 1, sizeof(double), fp );
	fread( (double *)&xUnit.y, 1, sizeof(double), fp );
	fread( (double *)&xUnit.z, 1, sizeof(double), fp );
	fread( (int *)&ySize, 1, sizeof(int), fp );
	fread( (double *)&yUnit.x, 1, sizeof(double), fp );
	fread( (double *)&yUnit.y, 1, sizeof(double), fp );
	fread( (double *)&yUnit.z, 1, sizeof(double), fp );
	fread( (int *)&zSize, 1, sizeof(int), fp );
	fread( (double *)&zUnit.x, 1, sizeof(double), fp );
	fread( (double *)&zUnit.y, 1, sizeof(double), fp );
	fread( (double *)&zUnit.z, 1, sizeof(double), fp );
	if( atom != NULL ) delete atom;
	atom = new CAtom[abs(IFlag)];
	for( int i=0; i< abs(IFlag); i++ ) atom[i].loadBin( fp );
	if( IFlag < 0 ){
		fread( (int *)&num1, 1, sizeof(int), fp );
		fread( (int *)&num2, 1, sizeof(int), fp );
	}
	if( d != NULL ) delete d;
	d = new float[xSize*ySize*zSize];
	fread( (float *)d, xSize*ySize*zSize, sizeof(float), fp );
	fclose( fp );
	return 0;
}

int CCube::saveBin( char *filename ){
FILE *fp;
	if( (fp = fopen( filename, "wb" )) == NULL ) return 1;
	fwrite( (int *)&IFlag, 1, sizeof(int), fp );
	fwrite( (double *)&origin.x, 1, sizeof(double), fp );
	fwrite( (double *)&origin.y, 1, sizeof(double), fp );
	fwrite( (double *)&origin.z, 1, sizeof(double), fp );
	fwrite( (int *)&xSize, 1, sizeof(int), fp );
	fwrite( (double *)&xUnit.x, 1, sizeof(double), fp );
	fwrite( (double *)&xUnit.y, 1, sizeof(double), fp );
	fwrite( (double *)&xUnit.z, 1, sizeof(double), fp );
	fwrite( (int *)&ySize, 1, sizeof(int), fp );
	fwrite( (double *)&yUnit.x, 1, sizeof(double), fp );
	fwrite( (double *)&yUnit.y, 1, sizeof(double), fp );
	fwrite( (double *)&yUnit.z, 1, sizeof(double), fp );
	fwrite( (int *)&zSize, 1, sizeof(int), fp );
	fwrite( (double *)&zUnit.x, 1, sizeof(double), fp );
	fwrite( (double *)&zUnit.y, 1, sizeof(double), fp );
	fwrite( (double *)&zUnit.z, 1, sizeof(double), fp );
	for( int i=0; i< abs(IFlag); i++ ) atom[i].saveBin( fp );
	if( IFlag < 0 ){
		fwrite( (int *)&num1, 1, sizeof(int), fp );
		fwrite( (int *)&num2, 1, sizeof(int), fp );
	}
	fwrite( (float *)d, xSize*ySize*zSize, sizeof(float), fp );
	fclose( fp );
	return 0;
}

int CCube::load( char *filename ){
	int i=0;
	while( filename[i] !='\0' ) i++;
	while( filename[i] !='.' ) i--;
	int result;
	if( filename[i+1] == 'c' ) result = loadCube( filename );
	else                       result = loadBin( filename );
	if( !result ) center = (xUnit * double(xSize-1) + yUnit * double(ySize-1) + zUnit * double(zSize-1)) * 0.5;
	return result;
}

int CCube::save( char *filename ){
	int i=0;
	while( filename[i] !='\0' ) i++;
	while( filename[i] !='.' ) i--;
	if( filename[i+1] == 'c' ) return saveCube( filename );
	else                       return saveBin( filename );
}