MCE - Marching Cube ELD SVN

#include "stdafx.h"
#include <math.h>  // matematic functions
#include <string.h>
#include <iostream>
#include <sstream>
#include <fstream>
#include "IO.h"
#include <iomanip>
#include <afxwin.h>
#include <afxinet.h>

//#include "VoxelMap.h"

IO::IO(void)
{
}
IO::~IO(void)
{
}
string IO::gen_random(int len) {
    //srand(time(0));
    static const char alphanum[] =
        "ABCDEFGHIJKLMNOPQRSTUVWXYZ";
    string s="";
    for (int i = 0; i < len; ++i) {
        s += alphanum[rand() % (sizeof(alphanum) - 1)];
    }
    return s;
}
bool IO::getFileFromURL(CString url, CString filename)
{
    CHttpFile *httpfile=NULL;
    try {		

        CInternetSession isession;
        httpfile = (CHttpFile*)isession.OpenURL(url, 1, INTERNET_FLAG_SECURE |
            INTERNET_FLAG_TRANSFER_BINARY | INTERNET_FLAG_RELOAD);
        if (httpfile)
        {
            DWORD dwStatusCode;
            httpfile->QueryInfoStatusCode(dwStatusCode);
            if(dwStatusCode == 200)
            {
                char buf[0x1000] = { 0 };
                DWORD read = 0;
                CFile file;
                if (file.Open(filename, CFile::modeCreate | CFile::modeWrite))
                    while(InternetReadFile(*httpfile, buf, sizeof(buf), &read) && read)
                        file.Write(buf, read);
            }

            httpfile->Close();
            delete httpfile;
        }
    } catch (CInternetException* pEx) {
      // catch any exceptions from WinINet      
		TCHAR szErr[1024];
		szErr[0] = '\0';        
        if(!pEx->GetErrorMessage(szErr, 1024))
			strcpy(szErr,"Some crazy unknown error");
		
		pEx->Delete();
		if(httpfile)
			delete httpfile;		
		return false;
	}
    return true;
}
bool IO::LoadDisplayStyleSettings(CString filename, DisplayStyleStruct &DS)
{
    CString strInputData;
	
	CFile Ff(filename, CFile::modeRead);
    if(!Ff) return false;
	CArchive ARar(&Ff,CArchive::load);

	while (ARar.ReadString(strInputData)) {
        vector<CString> values;
        strInputData.Replace("\t", " ");
        getValuesSeparatedByWhiteSpaces(strInputData, values);
        if(values.size()==2) {
            if(values[0].CompareNoCase("DisplayStyle")==0) {
                DS.DisplayStyle = (DisplayStyleMode) atoi(values[1]);
            } else if(values[0].CompareNoCase("BondsVisible")==0) {
                DS.bBondsVisible = atoi(values[1]);
            } else if(values[0].CompareNoCase("LineSize")==0) {
                DS.LineSize = atof(values[1]);
            } else if(values[0].CompareNoCase("BallSize")==0) {
                DS.BallSize = atof(values[1]);
            } else if(values[0].CompareNoCase("CylinderSize")==0) {
                DS.fCylinderSize = atof(values[1]);
            } else if(values[0].CompareNoCase("ScaledFactor")==0) {
                DS.ScaledFactor = atof(values[1]);
            } else if(values[0].CompareNoCase("BallResolution")==0) {
                DS.fBallResolution = atof(values[1]);
            } else if(values[0].CompareNoCase("CylinderResolution")==0) {
                DS.fCylinderResolution = atof(values[1]);
            } else if(values[0].CompareNoCase("ScaledBallResolution")==0) {
                DS.fScaledBallResolution = atof(values[1]);
            } else if(values[0].CompareNoCase("BondsColorStyle")==0) {
                DS.iBondsColorStyle = (BondsColorStyle) atoi(values[1]);
            } else if(values[0].CompareNoCase("Red")==0) {
                DS.BondsColor[0] = atof(values[1])/255;
            } else if(values[0].CompareNoCase("Green")==0) {
                DS.BondsColor[1] = atof(values[1])/255;
            } else if(values[0].CompareNoCase("blue")==0) {
                DS.BondsColor[2] = atof(values[1])/255;
            } else if(values[0].CompareNoCase("MapDisplayStyle")==0) {
                DS.MapDisplayStyle = (MapDisplayStyleMode) atoi(values[1]);
            } else if(values[0].CompareNoCase("MapLineSize")==0) {
                DS.fMapLineSize = atof(values[1]);
            } else if(values[0].CompareNoCase("MapCylinderSize")==0) {
                DS.fMapCylinderSize = atof(values[1]);
            } else if(values[0].CompareNoCase("MapCylinderResolution")==0) {
                DS.fMapCylinderResolution = atof(values[1]);
            } else if(values[0].CompareNoCase("UnitCellVisible")==0) {
                DS.bUnitCellVisible = atoi(values[1]);
            } else if(values[0].CompareNoCase("UnitCellLineSize")==0) {
                DS.fUnitCellLineSize = atof(values[1]);
            } else if(values[0].CompareNoCase("Organic")==0) {
                DS.bOrganic = atoi(values[1]);
            } else if(values[0].CompareNoCase("LoadMapSettings")==0) {
                DS.bLoadMapSettings = atoi(values[1]);
            } else if(values[0].CompareNoCase("MapsEdgesVisible")==0) {
                DS.bDrawMapEdges = atoi(values[1]);
            } else if(values[0].CompareNoCase("ProjectionMode")==0) {
                DS.PerspMode = (perspectiveMode) atoi(values[1]);
            } else if(values[0].CompareNoCase("PickingMode")==0) {
                DS.PickMode = (pickingMode) atoi(values[1]);
            } else if(values[0].CompareNoCase("MaxDistance")==0) {
                DS.maxdistance = atof(values[1]);
            }
        } else if(values.size()==4) {
            if(values[0].CompareNoCase("BackgroundColor")==0) {            
                DS.BgColor[0] = atof(values[1])/255;
                DS.BgColor[1] = atof(values[2])/255;
                DS.BgColor[2] = atof(values[3])/255;
            }            
        } else {
            //empty lines, etc...
        }
    };
    ARar.Close();
	Ff.Close();

    return true;
}
bool IO::SaveDisplayStyleSettings(CString filename, DisplayStyleStruct DS)
{
    char szBuffer[255];
    CFile f( filename, CFile::modeCreate | CFile::modeWrite);
    if(!f) return false;
	CArchive ar(&f,CArchive::store);
	
	ar.WriteString("		-- MCE Display style settings --		\n");
	ar.WriteString("		     -- generated by MCE --				\n\n");

	ar.WriteString("Molecule display style\n");
	sprintf(szBuffer,"DisplayStyle		%d\n", DS.DisplayStyle);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"BondsVisible		%d\n", DS.bBondsVisible);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"LineSize		%f\n", DS.LineSize);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"BallSize		%f\n", DS.BallSize);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"CylinderSize		%f\n", DS.fCylinderSize);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"ScaledFactor		%f\n", DS.ScaledFactor);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"BallResolution		%f\n", DS.fBallResolution);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"CylinderResolution	%f\n", DS.fCylinderResolution);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"ScaledBallResolution	%f\n", DS.fScaledBallResolution);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"BondsColorStyle		%d\n", DS.iBondsColorStyle);
	ar.WriteString(szBuffer);		
	sprintf(szBuffer,"Red			%f\n", DS.BondsColor[0]*255);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"Green			%f\n", DS.BondsColor[1]*255);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"blue			%f\n", DS.BondsColor[2]*255);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"MapDisplayStyle		%d\n", DS.MapDisplayStyle);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"MapLineSize		%f\n", DS.fMapLineSize);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"MapCylinderSize		%f\n", DS.fMapCylinderSize);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"MapCylinderResolution	%f\n", DS.fMapCylinderResolution);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"UnitCellVisible		%d\n", DS.bUnitCellVisible);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"UnitCellLineSize	%f\n", DS.fUnitCellLineSize);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"Organic	%d\n", DS.bOrganic);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"LoadMapSettings	%d\n", DS.bLoadMapSettings);
	ar.WriteString(szBuffer);
    sprintf(szBuffer,"MapsEdgesVisible	%d\n", DS.bDrawMapEdges);
	ar.WriteString(szBuffer);
    sprintf(szBuffer,"ProjectionMode	%d\n", DS.PerspMode);
	ar.WriteString(szBuffer);
    sprintf(szBuffer,"PickingMode	%d\n", DS.PickMode);
	ar.WriteString(szBuffer);
    sprintf(szBuffer,"MaxDistance	%f\n", DS.maxdistance);
	ar.WriteString(szBuffer);
    sprintf(szBuffer,"BackgroundColor	%f %f %f\n", DS.BgColor[0]*255, DS.BgColor[1]*255, DS.BgColor[2]*255);
	ar.WriteString(szBuffer);

	ar.Close();
	f.Close();
    return true;
}
bool IO::is_equal(float a, float b, float error) 
{
    if((a>(b-error)) && (a<(b+error))) return true;
    return false;
}

bool IO::isIdentityMatrix(float m[3][3])
{
    if(is_equal(m[0][0], 1.0) && is_equal(m[0][1], 0.0)  && is_equal(m[0][2], 0.0) &&
       is_equal(m[1][0], 0.0) && is_equal(m[1][1], 1.0)  && is_equal(m[1][2], 0.0) &&
       is_equal(m[2][0], 0.0) && is_equal(m[2][1], 0.0)  && is_equal(m[2][2], 1.0)) {
         return true;
    }
    return false;
}
bool IO::getRationalNumber(float n, CString &rational_number)
{
	rational_number = "";
	//get sign
	if((!is_equal(n, 0, 0.001)) && (n < 0)) {
    	rational_number = "-";
        n=-1*n;
    }
    
	//test 1,0,  1/2, 1/3, 2/3, 1/4, 2/4, 3/4, 1/6, 2/6, 3/6, 4/6, 5/6
	if(is_equal(n, 1, 0.001)) rational_number += "1";
	if(is_equal(n, 0, 0.001)) rational_number += "0";
	if(is_equal(n, 0.5, 0.001)) rational_number += "1/2";
	if(is_equal(n, 1.0/3.0, 0.001)) rational_number += "1/3";
	if(is_equal(n, 2.0/3.0, 0.001)) rational_number += "2/3";
	if(is_equal(n, 1.0/4.0, 0.001)) rational_number += "1/4";
	//if(is_equal(n, 2.0/4.0, 0.001)) rational_number += "2/4";
	if(is_equal(n, 3.0/4.0, 0.001)) rational_number += "3/4";
	if(is_equal(n, 1.0/6.0, 0.001)) rational_number += "1/6";
	//if(is_equal(n, 2.0/6.0, 0.001)) rational_number += "2/6";
	//if(is_equal(n, 3.0/6.0, 0.001)) rational_number += "3/6";
	//if(is_equal(n, 4.0/6.0, 0.001)) rational_number += "4/6";
	if(is_equal(n, 5.0/6.0, 0.001)) rational_number += "5/6";

	if((rational_number.GetLength()==0) || (rational_number.Compare("-")==0)) return false;
	else return true;
}
bool IO::getXYZRow(vector<CString> &row, CString &pos)
{
	if(row.size()!=4) return false;
	pos = "";
	if(row[0]=="1") pos += "x";
	if(row[0]=="-1") pos += "-x";
	if(row[1]=="1") { 
		if(pos.GetLength()!=0) pos += "+";
		pos += "y";
	}
	if(row[1]=="-1") pos += "-y";
	if(row[2]=="1") { 
		if(pos.GetLength()!=0) pos += "+";
		pos += "z";
	}
	if(row[2]=="-1") pos += "-z";
	
	if(row[3].Find("0")==-1) {
		if(row[3].Find("-")==-1) {
			pos += "+";
		}
		pos += row[3];
	}
	return true;
}
bool IO::ReadFouFile(CString filename, UnitCell *pUnit, VoxelMap *pMap)
{
	double PI = 3.1415926536;

	if(pUnit==0) return false;
	if(pMap==0) return false;

	CString strInputData;

	// support variable
	int nLineNumber=0;
	int i, j;
	int x,y,z;
	
	CFile f;
	if(!f.Open(filename, CFile::modeRead)) {
		AfxMessageBox("File: " + filename + " does not exists."); 
		return false;
	}
	CArchive ar(&f,CArchive::load);

	
	ar.ReadString(strInputData); //INFO DOWN, ....?
	ar.ReadString(strInputData); //TRAN ?
	 
    float m[3][3], mi[3][3];
	for (i=0;i<3;i++)
	{
		for(j=0;j<3;j++)
		{
			ar.ReadString(strInputData);
			m[i][j]=(float) atof(strInputData);
			TRACE( "%d %d %f \n",i,j,m[i][j]);
		}
	}
    
    for (i=0;i<3;i++)
	{
		ar.ReadString(strInputData);
		pUnit->Translation_Matrix[i]=(float) atof(strInputData);
	}

	ar.ReadString(strInputData); // CELL
	
	ar.ReadString(strInputData);
	pUnit->A = (float) atof(strInputData);
	ar.ReadString(strInputData);
	pUnit->B = (float) atof(strInputData);
	ar.ReadString(strInputData);
	pUnit->C = (float) atof(strInputData);
	ar.ReadString(strInputData);
	pUnit->Alpha = (float) (atof(strInputData)/PI*180);
	ar.ReadString(strInputData);
	pUnit->Beta = (float) (atof(strInputData)/PI*180);
	ar.ReadString(strInputData);
	pUnit->Gama = (float) (atof(strInputData)/PI*180);

    //if there is just identity matrix, voxels are in the whole unit cell
    if(isIdentityMatrix(m)) {        
        pMap->mType = GSAS_GRD_FILE; // Map type identification
        pUnit->SetSpaceGroup("P1");

        float m1[3][3], m2[3][3];
        pUnit->getFTC_Matrix(m1);
        pUnit->getCTF_Matrix(m2);
	    for(i=0;i<3;i++) {
		    for(j=0;j<3;j++){
                pMap->MatrixFtC[i][j] = m1[i][j];
                pMap->MatrixCtF[i][j] = m2[i][j];
		    }
        }
    } else {
        pMap->mType = CRYSTALS_FOU_FILE; // Map type identification
        UnitCell::InvertMatrix(m, mi);
        pUnit->setFtC_Matrix(m); 
    }

	ar.ReadString(strInputData); // L14

    /*
    	for (i=0;i<3;i++)
	{
		pMap->AxisParam[i].Start=0;		
		pMap->AxisParam[i].End=1;
		pMap->AxisParam[i].Division=stepsx;
		
	}
	pMap->AxisParam[0].Step=a/((float)stepsx-1);
	pMap->AxisParam[1].Step=b/((float)stepsy-1);
	pMap->AxisParam[2].Step=c/((float)stepsz-1);
   	pMap->AxisParam[0].Dimension=a;
	pMap->AxisParam[1].Dimension=b;
	pMap->AxisParam[2].Dimension=c;


	pMap->Lines=stepsx;
	pMap->Columns=stepsy;
	pMap->Sections=stepsz;
    */
	
	for (i=0;i<3;i++)
	{
		ar.ReadString(strInputData);
		pMap->AxisParam[i].Start=(float) atof(strInputData);
        if(pMap->mType == GSAS_GRD_FILE) {
            pMap->AxisParam[i].Start = 0;
        }

		ar.ReadString(strInputData);
		pMap->AxisParam[i].Step=(float) atof(strInputData);

		ar.ReadString(strInputData);
		pMap->AxisParam[i].End=(float) atof(strInputData);
        if(pMap->mType == GSAS_GRD_FILE) {
            pMap->AxisParam[i].End = 1;
        }

		ar.ReadString(strInputData);
		pMap->AxisParam[i].Division=(float) atof(strInputData);
        
	}
    pMap->AxisParam[0].Dimension = pUnit->A;
    pMap->AxisParam[1].Dimension = pUnit->B;
    pMap->AxisParam[2].Dimension = pUnit->C;

	ar.ReadString(strInputData); // SIZE

	ar.ReadString(strInputData);
	pMap->Lines=atoi(strInputData);

	ar.ReadString(strInputData);
	pMap->Columns=atoi(strInputData);

	ar.ReadString(strInputData);
	pMap->Sections=atoi(strInputData);

	if ((pMap->Lines>MAX_OPTIMAL_VOXEL_DIEMENSION)|
		(pMap->Columns>MAX_OPTIMAL_VOXEL_DIEMENSION)|
		(pMap->Sections>MAX_OPTIMAL_VOXEL_DIEMENSION)){
		AfxMessageBox("Too large ELD field(in one direction bigger than 400)!\n Program will use more memory...\n");
	}

	pMap->MemoryForVoxel();
    
	
	// !!! THis part of code is not memory-optimal !!!

	for (z=0;z<pMap->Sections;z++)
	{
		ar.ReadString(strInputData); // BLOCK
		ar.ReadString(strInputData); // N

		for(x=0;x<pMap->Lines;x++)
		{
			for(y=0;y<pMap->Columns;y++)
			{
				ar.ReadString(strInputData);
				pMap->Voxel[x][y][z]=(float) atof(strInputData);
			}
		}
	}
	
	TRACE("Finished reading Voxel Map");

	// mlecular data input

	ar.ReadString(strInputData); //LIST5

	ar.ReadString(strInputData);
	int nbAtoms = atoi(strInputData);

	//if(nbAtoms>IO_MAXATOMS) {
		//AfxMessageBox("Error: There are more than 2000 atoms in input file.");
		//return false;
	//}

	if (nbAtoms>0)
	{
		ar.ReadString(strInputData);
		int atom_record_lenght=atoi(strInputData);
        switch (atom_record_lenght) {
            case 14: {
                for (i=0;i<nbAtoms;i++) {
			        ar.ReadString(strInputData); // atom name
			        CString AtName = strInputData.Trim();
			        ar.ReadString(strInputData);
			        int AtNumber = atoi(strInputData);
                    CString AtLabel;
                    AtLabel.Format("%s(%d)", AtName, AtNumber); 
			        ar.ReadString(strInputData);
			        float occ =(float) atof(strInputData);
			        ar.ReadString(strInputData); // reserved, dummy now
			        ar.ReadString(strInputData);
			        float Atx =(float) atof(strInputData);
			        ar.ReadString(strInputData);
			        float Aty=(float) atof(strInputData);
			        ar.ReadString(strInputData);
			        float Atz=(float) atof(strInputData);            
                    AtName.Trim();
                    if(pMap->mType == GSAS_GRD_FILE) {
                        pUnit->CalcCartToFract(Atx, Aty, Atz);
                        Atx+=0.5;
                        Aty+=0.5;
                        Atz+=0.5;
                    }
			        pUnit->AddNewAtomInFract(AtName, AtLabel, Atx, Aty, Atz); 
			        ar.ReadString(strInputData);
			        float uiso =(float) atof(strInputData);
			        for(j=0;j<6;j++) {
				        ar.ReadString(strInputData);
				        float uaniso=(float) atof(strInputData);
			        }			
                }
                pUnit->SetSpaceGroupSpecial("P1");
             }
             break;
            case 18: {
                for (i=0;i<nbAtoms;i++) {
			        ar.ReadString(strInputData); // atom name
			        CString AtName = strInputData.Trim();
			        ar.ReadString(strInputData);
			        int AtNumber = atoi(strInputData);
                    CString AtLabel;
                    AtLabel.Format("%s(%d)", AtName, AtNumber); 
			        ar.ReadString(strInputData);
			        float occ =(float) atof(strInputData);
			        ar.ReadString(strInputData); // reserved, dummy now
			        ar.ReadString(strInputData);
			        float Atx =(float) atof(strInputData);
			        ar.ReadString(strInputData);
			        float Aty=(float) atof(strInputData);
			        ar.ReadString(strInputData);
			        float Atz=(float) atof(strInputData);            
                    AtName.Trim();
                    if(pMap->mType == GSAS_GRD_FILE) {
                        pUnit->CalcCartToFract(Atx, Aty, Atz);
                        Atx+=0.5;
                        Aty+=0.5;
                        Atz+=0.5;
                    }
			        pUnit->AddNewAtomInFract(AtName, AtLabel, Atx, Aty, Atz);
			        ar.ReadString(strInputData);
			        float uiso =(float) atof(strInputData);
			        for(j=0;j<6;j++) {
				        ar.ReadString(strInputData);
				        float uaniso=(float) atof(strInputData);
			        }
			        ar.ReadString(strInputData);
                    ar.ReadString(strInputData);
                    ar.ReadString(strInputData);
                    ar.ReadString(strInputData);
                }
                pUnit->SetSpaceGroupSpecial("P1");
            }
            break;
        case 7:{
                for (i=0;i<nbAtoms;i++) {
                    ar.ReadString(strInputData); // atom name
			        CString AtName = strInputData.Trim();
			        ar.ReadString(strInputData);
			        int AtNumber = atoi(strInputData);
                    CString AtLabel;
                    AtLabel.Format("%s(%d)", AtName, AtNumber); 
			        ar.ReadString(strInputData);			        
			        float Atx =(float) atof(strInputData);
			        ar.ReadString(strInputData);
			        float Aty=(float) atof(strInputData);
			        ar.ReadString(strInputData);
			        float Atz=(float) atof(strInputData);            
                    AtName.Trim();
                    if(pMap->mType == GSAS_GRD_FILE) {
                        pUnit->CalcCartToFract(Atx, Aty, Atz);
                        Atx+=0.5;
                        Aty+=0.5;
                        Atz+=0.5;
                    }
			        pUnit->AddNewAtomInFract(AtName, AtLabel, Atx, Aty, Atz);
			        ar.ReadString(strInputData);//dummy
			        ar.ReadString(strInputData);//dummy                     
                } 
                //pUnit->SetSpaceGroup("P1");

                /*float m[3][3], mi[3][3];
	            /for (i=0;i<3;i++)
	            {
		            for(j=0;j<3;j++)
		            {			            
			            m[i][j]=0;                        
		            }
	            }
                m[0][0] = pUnit->A;
                m[1][1] = pUnit->B;
                m[2][2] = pUnit->C;
                UnitCell::InvertMatrix(m, mi);
                pUnit->setFtC_Matrix(m);  
                */

                pUnit->SetSpaceGroupSpecial("P1");
               }
            break;
        default: {
             for (i=0;i<nbAtoms;i++) {
			        ar.ReadString(strInputData); // atom name
			        CString AtName = strInputData.Trim();
			        ar.ReadString(strInputData);
			        int AtNumber = atoi(strInputData);
                    CString AtLabel;
                    AtLabel.Format("%s(%d)", AtName, AtNumber);
			        ar.ReadString(strInputData);
			        float occ =(float) atof(strInputData);
			        ar.ReadString(strInputData); // reserved, dummy now
			        ar.ReadString(strInputData);
			        float Atx =(float) atof(strInputData);
			        ar.ReadString(strInputData);
			        float Aty=(float) atof(strInputData);
			        ar.ReadString(strInputData);
			        float Atz=(float) atof(strInputData);            
                    AtName.Trim();
                    if(pMap->mType == GSAS_GRD_FILE) {
                        pUnit->CalcCartToFract(Atx, Aty, Atz);
                        Atx+=0.5;
                        Aty+=0.5;
                        Atz+=0.5;
                    }
			        pUnit->AddNewAtomInFract(AtName, AtLabel, Atx, Aty, Atz);
			        ar.ReadString(strInputData);
			        float uiso =(float) atof(strInputData);
			        for(j=0;j<6;j++) {
				        ar.ReadString(strInputData);
				        float uaniso=(float) atof(strInputData);
			        }			
			        // dumy data from future fou versions
			        if (atom_record_lenght>14)
			        {
				        for(j=15;j<=atom_record_lenght;j++)
				        {
					        ar.ReadString(strInputData);
				        }
			        }
                }
             pUnit->SetSpaceGroupSpecial("P1");
            }
            break;
        }	
    } // endif atom presents
	TRACE("Finished reading Molecule data");

	ar.Close();
	f.Close();
    

	return true;
}


bool IO::ReadGSASFobFile(CString filename, UnitCell *pUnit, VoxelMap *pMap)
{
	FILE *fouin;
	int dummy_int, nRead, i, interactive = 1;
	long int record_size, records_per_section;
	FourierHeader hdr;
	FourierData *data;
	PhaseData phase_data;

	if(pUnit==0) return false;
	if(pMap==0) return false;


	pMap->mType = GSAS_FOB_FILE;
	if ((fouin = fopen(filename, "rb")) != NULL)
	{	
		// Read in header information.  
		fread(hdr.phase, sizeof(char), 66, fouin); hdr.phase[66] = '\0';
		fread(hdr.title, sizeof(char), 66, fouin); hdr.title[66] = '\0';
		fread(hdr.maptyp, sizeof(char), 4, fouin); hdr.maptyp[4] = '\0';

		do {
			nRead = fread(&dummy_int, sizeof(int), 1, fouin);
		} while (dummy_int == 0 && nRead == 1);


		record_size = ftell(fouin) - sizeof(int);
		hdr.nxyzi[0] = dummy_int;
		fread(&hdr.nxyzi[1], sizeof(int), 2, fouin);
		fread(hdr.nxyzo, sizeof(int), 3, fouin);
		fread(hdr.nxyzt, sizeof(int), 3, fouin);
		fread(&hdr.msect, sizeof(int), 1, fouin);
		fread(&hdr.nrho, sizeof(int), 1, fouin);
		fread(&hdr.rmax, sizeof(float), 1, fouin);
		fread(&hdr.rmin, sizeof(float), 1, fouin);
		fread(&hdr.srho, sizeof(float), 1, fouin);
		if ((hdr.nrho + 2)*4 > record_size) {
			records_per_section = ((hdr.nrho + 2)*4)/record_size + 1;
		} else {
			records_per_section = 1;
		}

		if (hdr.msect > 3 || hdr.msect < 1) {
			MessageBox(NULL, "Bad msect value, hdr.msect. Should be 1, 2 or 3.\n Are you sure that the map file was\n created on this machine?", "WinMain", MB_ICONEXCLAMATION); 
			return FALSE;
		}

		data = (FourierData*) malloc(hdr.nxyzt[2]*sizeof(FourierData));
		if (data == NULL) {
			MessageBox(NULL, "Error allocating memory for Fourier data.\n", "WinMain", MB_ICONEXCLAMATION); 
			//fprintf(stderr, "Error allocating memory for Fourier data.\n");
			fclose(fouin);
			exit(1);
		}
		for (i = 0; i < hdr.nxyzt[2]; i++) {
			data[i].rho = (float *) malloc(hdr.nrho*sizeof(float));
			if (data[i].rho == NULL) {
				MessageBox(NULL, "Error allocating memory for Fourier sections.\n", "WinMain", MB_ICONEXCLAMATION); 
				free(data);
				fclose(fouin);
				exit(1);
				}
		}

		for (i = 0; i < hdr.nxyzt[2]; i++) {
			// Jump to start of section i. 
			fseek(fouin, record_size*(2 + records_per_section*i), SEEK_SET);
			fread(data[i].rho, sizeof(float), hdr.nrho, fouin);
			fread(&data[i].romx, sizeof(float), 1, fouin);
			fread(&data[i].romn, sizeof(float), 1, fouin);    
		}
		fclose(fouin);

		ReadGSASFobFile_parse_expfile(filename, &phase_data, interactive);

		ReadGSASFobFile_output_crystals(&phase_data, &hdr, data, interactive, pUnit, pMap);
		
		pUnit->SetSpaceGroupSpecial("P1");
		//pUnit->InvertMatrix(pUnit->FtC_Matrix, pUnit->CtF_Matrix);

		for (i = 0; i < hdr.nxyzt[2]; i++) free(data[i].rho);
		free(data);
	} else {	
		AfxMessageBox("File: " + filename + " does not exists."); 
		return false;
	}
	return true;
}
bool IO::Read3DFile(CString filename, UnitCell *pUnit, VoxelMap *pMap)
{
    if(pUnit==0) return false;
	if(pMap==0) return false;

	CString strInputData;
	int nLineNumber=0;
	int i, j, k;
	int x,y,z;
	CFile f;
	if(!f.Open(filename, CFile::modeRead)) {
		AfxMessageBox("File: " + filename + " does not exists."); 
		return false;
	}
	CArchive ar(&f,CArchive::load);
    pMap->mType = THREE_D;
    /*
   11.8467     11.8467      7.5032     90.00000    90.00000    90.00000
   116     116      76
     0       0       0         0.0022200
     1       0       0         0.0022400
     2       0       0         0.0023100

    */
    //first line
    float a,b,c,al,be,ga;
    ar.ReadString(strInputData); // unit cell parameters
	sscanf(strInputData,"%f %f %f %f %f %f",&a,&b,&c,&al,&be,&ga);
	pUnit->A = a;
	pUnit->B = b;
	pUnit->C = c;
	pUnit->Alpha = al;
	pUnit->Beta = be;
	pUnit->Gama = ga;

    //second line
    int N_x, N_y, N_z;
	ar.ReadString(strInputData); // steps in lattice directions
	sscanf(strInputData,"%d %d %d",&N_x, &N_y, &N_z);

	// let simulate start at 0 and 0.1 strps    
    //The density is the value at the center of the voxel / pixel!!!
	for (i=0;i<3;i++)
	{
		pMap->AxisParam[i].Start=0;
		pMap->AxisParam[i].Step=0.1f;
		pMap->AxisParam[i].End=1;
		pMap->AxisParam[i].Division=0.1f;
	}
    pMap->AxisParam[0].Dimension=a;
    pMap->AxisParam[1].Dimension=b;
    pMap->AxisParam[2].Dimension=c;

	pMap->Lines=N_x;
	pMap->Columns=N_y;
	pMap->Sections=N_z;

	if ((pMap->Lines>MAX_OPTIMAL_VOXEL_DIEMENSION)|
		(pMap->Columns>MAX_OPTIMAL_VOXEL_DIEMENSION)|
		(pMap->Sections>MAX_OPTIMAL_VOXEL_DIEMENSION)){
		AfxMessageBox("Too large ELD field(in one direction bigger than 400)!\n Program will use more memory...\n");
	}

	pMap->MemoryForVoxel();
    float rho;
    while(ar.ReadString(strInputData)) {
        sscanf(strInputData,"%d %d %d %f",&i, &j, &k, &rho);
        if(i<0 || j<0 || k<0) {
            AfxMessageBox("Error in reading file i, j or k index is < 0\n");
            return false;
        }
        if((i<N_x) && (j<N_y) && (k<N_z)) {
            pMap->Voxel[i][j][k] = rho;
        } else {
            AfxMessageBox("Error in reading file i, j or k index is bigger than dimension\n");
            return false;
        }
    };

    /*
	for(x=0;x<pMap->Lines;x++)
	{
		for(y=0;y<pMap->Columns;y++)
		{
			for (z=0;z<pMap->Sections;z++)
			{
				pMap->Voxel[x][y][z] = GetNumberFromFile(&ar);
			}
		}
	}
    */

	TRACE("Finished reading Voxel Map");

	ar.Close();
	f.Close();	


    pUnit->SetSpaceGroup("P1");
    float m1[3][3], m2[3][3];
    pUnit->getFTC_Matrix(m1);
    pUnit->getCTF_Matrix(m2);
	for(i=0;i<3;i++) {
		for(j=0;j<3;j++){
            pMap->MatrixFtC[i][j] = m1[i][j];
            pMap->MatrixCtF[i][j] = m2[i][j];
		}
    }

	return true;
}
bool IO::ReadXplorFile(CString filename, UnitCell *pUnit, VoxelMap *pMap)
{
    if(pUnit==0) return false;
	if(pMap==0) return false;

	CString strInputData;
	CFile f;
	if(!f.Open(filename, CFile::modeRead)) {
		AfxMessageBox("File: " + filename + " does not exists."); 
		return false;
	}
	CArchive ar(&f,CArchive::load);

	ar.ReadString(strInputData); //read empty line

    ar.ReadString(strInputData); // line 1 of the title
    int title_line;
    sscanf(strInputData,"%d",&title_line);
    for(int i=0;i<title_line;i++) {
        ar.ReadString(strInputData); // title
    }
	
    int x_div, x_start, x_end, y_div, y_start, y_end, z_div, z_start, z_end;
    ar.ReadString(strInputData); //pixel division: number of pixels along x, first pixel, last pixel; ditto for y and z    
    sscanf(strInputData,"%d %d %d %d %d %d %d %d %d", &x_div, &x_start, &x_end, &y_div, &y_start, &y_end, &z_div, &z_start, &z_end);

    ar.ReadString(strInputData); //cell parameters
	float a,b,c,al,be,ga;
	sscanf(strInputData,"%f %f %f %f %f %f",&a,&b,&c,&al,&be,&ga);
	pUnit->A = a;
	pUnit->B = b;
	pUnit->C = c;
	pUnit->Alpha = al;
	pUnit->Beta = be;
	pUnit->Gama = ga;

    CString order;
    ar.ReadString(order); // ZYX string defining the order of axes from the most slowly to the most quickly varying

	for (int i=0;i<3;i++)
	{
		pMap->AxisParam[i].Start=0;
		pMap->AxisParam[i].Step=0.1f;
		pMap->AxisParam[i].End=1;
		pMap->AxisParam[i].Division=0.1f;
	}
    pMap->AxisParam[0].Dimension=a;
    pMap->AxisParam[1].Dimension=b;
    pMap->AxisParam[2].Dimension=c;

	pMap->Lines=x_div;
	pMap->Columns=y_div;
	pMap->Sections=z_div;

	if ((pMap->Lines>MAX_OPTIMAL_VOXEL_DIEMENSION)|
		(pMap->Columns>MAX_OPTIMAL_VOXEL_DIEMENSION)|
		(pMap->Sections>MAX_OPTIMAL_VOXEL_DIEMENSION)){
		AfxMessageBox("Note: Too large ELD field(in one direction bigger than 400)!\n Program will use more memory...\n");
	}

	pMap->MemoryForVoxel();

	for(int x=0;x<pMap->Lines;x++) {
		for(int y=0;y<pMap->Columns;y++) {
			for (int z=0;z<pMap->Sections;z++) {
				pMap->Voxel[x][y][z] = 0; //GetNumberFromFile(&ar);
			}
		}
	}
    if(order.Compare("ZYX")==0) {
        for(int z=z_start;z<=z_end;z++) {
            int test;
            ar.ReadString(strInputData);
            sscanf(strInputData,"%d",&test);
            CString line="";
		    for(int y=y_start;y<=y_end;y++) {
			    for (int x=x_start;x<=x_end;x++) {
				    if(!GetNumberFromString(line, 12, pMap->Voxel[x][y][z])){
                        ar.ReadString(line);
                        GetNumberFromString(line, 12, pMap->Voxel[x][y][z]);
                    }
			    }
		    }            
	    }
    } else if(order.Compare("ZXY")==0) {
        for(int z=z_start;z<=z_end;z++) {
            int test;
            ar.ReadString(strInputData);
            sscanf(strInputData,"%d",&test);
            CString line="";
		    for(int x=x_start;x<=x_end;x++) {
			    for (int y=y_start;y<=y_end;y++) {
				    if(!GetNumberFromString(line, 12, pMap->Voxel[x][y][z])){
                        ar.ReadString(line);
                        GetNumberFromString(line, 12, pMap->Voxel[x][y][z]);
                    }
			    }
		    }            
	    }
    } else if(order.Compare("XZY")==0) {
        for(int x=x_start;x<=x_end;x++) {
            int test;
            ar.ReadString(strInputData);
            sscanf(strInputData,"%d",&test);
            CString line="";
		    for(int z=z_start;z<=z_end;z++) {
			    for (int y=y_start;y<=y_end;y++) {
				    if(!GetNumberFromString(line, 12, pMap->Voxel[x][y][z])){
                        ar.ReadString(line);
                        GetNumberFromString(line, 12, pMap->Voxel[x][y][z]);
                    }
			    }
		    }            
	    }
    } else if(order.Compare("XYZ")==0) {
        for(int x=x_start;x<=x_end;x++) {
            int test;
            ar.ReadString(strInputData);
            sscanf(strInputData,"%d",&test);
            CString line="";
		    for(int y=y_start;y<=y_end;y++) {
			    for (int z=z_start;z<=z_end;z++) {
				    if(!GetNumberFromString(line, 12, pMap->Voxel[x][y][z])){
                        ar.ReadString(line);
                        GetNumberFromString(line, 12, pMap->Voxel[x][y][z]);
                    }
			    }
		    }            
	    }
    } else {
        AfxMessageBox("Error while reading the file: '" + order + "' order is not defined.\n Map was not loaded."); 
    }

    //xplor reads map from 0 to (1-step) of the unit cell...
    //complete map in the whole unit cell from 0 to 1.
    pMap->complete_map_in_direction(nX);
    pMap->complete_map_in_direction(nY);
    pMap->complete_map_in_direction(nZ);

	TRACE("Finished reading Voxel Map");
  
	ar.Close();
	f.Close();

	pUnit->SetSpaceGroup("P1");
    float m1[3][3], m2[3][3];
    pUnit->getFTC_Matrix(m1);
    pUnit->getCTF_Matrix(m2);
	for(int i=0;i<3;i++)
		for(int j=0;j<3;j++){
            pMap->MatrixFtC[i][j] = m1[i][j];
            pMap->MatrixCtF[i][j] = m2[i][j];
		}

    pMap->isPeriodic = true;
	return true;
}
bool IO::ReadGRDFile(CString filename, UnitCell *pUnit, VoxelMap *pMap)
{
	if(pUnit==0) return false;
	if(pMap==0) return false;

	CString strInputData;
	CFile f;
	if(!f.Open(filename, CFile::modeRead)) {
		AfxMessageBox("Can't open file: " + filename + "."); 
		return false;
	}
	CArchive ar(&f,CArchive::load);

	ar.ReadString(strInputData); 
	
	if(strInputData == "2DGRDFIL  0"){
		pMap->mType = XD_GRD2D_FILE;
        pMap->m_2D_map.setVisibility(true);
	}
	if (strInputData == "3DGRDFIL  0"){
		pMap->mType = XD_GRD3D_FILE;        
	}
	ar.Close();
	f.Close();

	if((pMap->mType==XD_GRD2D_FILE)|(pMap->mType==XD_GRD3D_FILE)) {
		ReadXDGRDFile(filename, pUnit, pMap);        
    }
	else {
		pMap->mType = GSAS_GRD_FILE;
        pMap->isPeriodic = true;
		ReadNormalGRDFile(filename, pUnit, pMap);
	}

	return true;
}
bool IO::ImportFragmentFromCif(CString filename, Molecule *mol)
{   
    if(mol==NULL) return false;

    CifReader cif;

    if(!cif.readFile(filename.GetString())) {
        AfxMessageBox("Error during file reading (" + filename + ")!"); 
        return false;
    }

    cif_crystal crystal = cif.getCrystal();
    std::vector<cif_atom> atoms = cif.getAtoms();

    UnitCell tmpUnit;
    tmpUnit.A = crystal.a;
    tmpUnit.B = crystal.b;
    tmpUnit.C = crystal.c;
    tmpUnit.Alpha = crystal.al;
    tmpUnit.Beta = crystal.be;
    tmpUnit.Gama = crystal.ga;
    tmpUnit.SetSpaceGroup(crystal.sg.c_str());
    if(crystal.equiv_pos.size()!=0) {
        tmpUnit.listMatrix = crystal.equiv_pos_matrixes;        
    }

    for(int i=0;i<atoms.size();i++) {
         tmpUnit.AddNewAtomInFract(atoms[i].type.c_str(), atoms[i].name.c_str(), atoms[i].x, atoms[i].y, atoms[i].z);
    }

    std::vector<Atom> fragment = tmpUnit.getAssym_part_atoms_cart();

    for(int i=0;i<fragment.size();i++) {
        mol->AddNewAtom(fragment[i]);
    }
    
    /*
    for(int i=0;i<fragment.size();i++) {
        pUnit->CalcCartToFract(&fragment[i]);   
    }

    //move to O, O, O position
    Molecule mol;
    for(int i=0;i<fragment.size();i++) {
        mol.AddNewAtom(fragment[i]);
    }
    float center[3];
    mol.GetCenter(center[0], center[1], center[2]);
    center[0] = -1*center[0];
    center[1] = -1*center[1];
    center[2] = -1*center[2];

    mol.Move(center);

    int nb = pUnit->getAssym_part_molecule_fract().size();
    for(int i=0;i<mol.GetAtomNb();i++) {
        pUnit->AddNewAtomInFract(mol.GetAtom(i), nb, false, true);
    }    
    */
    return true;
}
std::string IO::ltrim(std::string s)
{
    int i;
    for(i=0;i<s.length();i++) {
        if (!isspace(s.at(i))) {
            break;
        } 
    }
    return s.substr(i, s.length());
}
std::string IO::rtrim(std::string s)
{
    int i;
    for(i=s.length()-1;i>=0;i--) {
        if (!isspace(s.at(i))) {
            break;
        } 
    }
    return s.substr(0, i+1);
}
bool IO::getValuesSeparatedByWhiteSpaces(CString list, vector<CString> &res)
{	
    list.Trim();
	if(list.GetLength()==0) return true;
	list+=" ";
	do {
		int pos = list.Find(" ");
		res.push_back(list.Mid(0, pos));
		list = list.Mid(pos, list.GetLength()-pos);
        list.TrimLeft();
	} while(list.GetLength()>0);

	return true; 
}
Molecule IO::getMoleculeFromFragmentDBEntry(FragmentDBStruct entry)
{
    UnitCell unit(entry.a, entry.b, entry.c, entry.al, entry.be, entry.ga);
    unit.GenerateFractToOrtMatrix();
    Molecule mol;
    for(int i=0;i<entry.atoms.size();i++) {
        Atom tmpat = entry.atoms[i];
        unit.CalcFractToCart(&tmpat);
        mol.AddNewAtom(tmpat);
    }   
    float c[3] = {0,0,0};
    //mol.GetCenter(c[0], c[1], c[2]);
    mol.MoveToPoint(c);
    return mol;
}
bool IO::SaveFragmentDBFile(CString filename, vector<FragmentDBStruct> &FDB)
{
    char szBuffer[250];
    //saving the old file - make backup
    if(PathFileExistsA(filename)) {
        if(PathFileExistsA(filename+".bak")) {
            DeleteFileA(filename+".bak");
        }
        MoveFile(filename, filename+".bak");    
    }
    

	CFile f( filename, CFile::modeCreate | CFile::modeWrite);
	CArchive ar(&f,CArchive::store);
	
    ar.WriteString("# Fragment database of DSR.py\n");
    ar.WriteString("#\n");
    ar.WriteString("# Some Fragments are from geometry optimizations with Gaussian 03:\n");
    ar.WriteString("#  Gaussian 03, Revision B.04,\n");
    ar.WriteString("#  M. J. Frisch, et. al. Gaussian, Inc., Pittsburgh PA, 2003.\n");
    ar.WriteString("#\n");
    ar.WriteString("# Some with TURBOMOLE V6.0 by K.Eichkorn, O.Treutler, H.Oehm,\n");
    ar.WriteString("#  M.Haeser and R.Ahlrichs (Chemical Physics Letters 242 (1995) 652-660)\n");
    ar.WriteString("#  parallel version: M.v.Arnim  &  R.Ahlrichs quantum chemistry group\n");
    ar.WriteString("#  university  karlsruhe germany\n");
    ar.WriteString("#\n");
    ar.WriteString("# Some are from Ilia A. Guzeis Idealized Molecular Geometry Library:\n");
    ar.WriteString("#  Guzei, Ilia A. 2010. Idealized Molecular Geometry Library.\n");
    ar.WriteString("#  University of Wisconsin-Madison. Madison, WI, USA\n");
    ar.WriteString("#  http://xray.chem.wisc.edu/Projects/IdealizedMolecularGeometry.html \n");
    ar.WriteString("#\n");
    ar.WriteString("# Some are from the Cambridge structural Database:\n");
    ar.WriteString("# Groom, C. R. & Allen, F. H. (2014). Angew. Chem. Int. Ed. 53, 662-671.\n");
    ar.WriteString("# https://www.ccdc.cam.ac.uk/\n");
    ar.WriteString("#\n");
    ar.WriteString("# Some are imported from the Grade Web Server http://grade.globalphasing.org\n");
    ar.WriteString("#\n");
    ar.WriteString("# Other fragments are from in-house measurements.\n");
    ar.WriteString("#\n");
    ar.WriteString("########################################################################################\n");
    ar.WriteString("#   You can insert your own fragments to \"dsr_user_db.txt\". Just follow this syntax:\n");
    ar.WriteString("########################################################################################\n");
    ar.WriteString("# Database syntax:\n");
    ar.WriteString("# <name>                           <- Start tag.\n");
    ar.WriteString("# RESI class                       <- Required, defines the residue name of db entry.\n");
    ar.WriteString("# RESTRAINTS                       <- Any restraints and comments following the SHELXL syntax.\n");
    ar.WriteString("#                                     You must enter at least one restraint! e.g. RIGU C1 > C7\n");
    ar.WriteString("# FRAG 17 a b c alpha beta gamma   <- FRAG card with afix number and cell parameters.\n");
    ar.WriteString("#      Atom  number  coordinates   <- One isotropic atom per line. \n");
    ar.WriteString("# e.g.  O1     1       x  y  z     <- Either the atom type is recognized by the atom name for \n");
    ar.WriteString("#                                     positive numbers.\n");
    ar.WriteString("# e.g.  O1    -8       x  y  z     <- Or the atom type is defined by the negative atomic     \n");
    ar.WriteString("#                                     number.\n");
    ar.WriteString("# </name>                          <- End tag. Same as start tag but with /\n");
    ar.WriteString("#\n");
    ar.WriteString("# - Anything not being an atom after FRAG is ignored.\n");
    ar.WriteString("# - Fragment names CF3, CF6 and CF9 are reserved by DSR. Do not attempt to use them in any\n");
    ar.WriteString("#   database entry.\n");
    ar.WriteString("# - Only lines beginning with valid SHELXL instructions are allowed in the header.\n");
    ar.WriteString("# - Anything behind the 5th column in the atom list is ignored.\n");
    ar.WriteString("# - Long lines can be wrapped with an equal sign (=) and a space character in the next line\n");
    ar.WriteString("#   like in SHELXL, but they can also be of any length. All lines will be wrapped to fit\n");
    ar.WriteString("#   in the SHELXL file automatically.\n");
    ar.WriteString("########################################################################################\n\n");

    for(int i=0;i<FDB.size();i++) {
        ar.WriteString("<"+FDB[i].tag_name+">\n");                
        //ar.WriteString("rem Name: "+FDB[i].alt_names+"\n"); - it is part of the rem...
        ar.WriteString("REM Name: "+FDB[i].alt_names.Trim()+"\n");
        ar.WriteString("REM Src: "+FDB[i].src.Trim()+"\n");
        if(FDB[i].rem.Trim().GetLength()!=0) {
            ar.WriteString(FDB[i].rem.Trim()+"\n");
        }
        if(FDB[i].resi.Trim().GetLength()!=0) {
            ar.WriteString(FDB[i].resi.Trim()+"\n");
        }
        FDB[i].shelx_restr.Replace("\r\n", "\n");
        if(FDB[i].shelx_restr.Trim().GetLength()!=0) {
            ar.WriteString(FDB[i].shelx_restr.Trim()+"\n");
        }
        CString lines = FDB[i].jana_restr;
        CString line;
        while(IO::ReadLineFromString(lines, line)) {
            ar.WriteString("JANA_RESTR "+line+"\n");
        }
        lines = FDB[i].Crystals_restr;
        while(IO::ReadLineFromString(lines, line)) {
            ar.WriteString("CRYSTALS_RESTR "+line+"\n");
        }
        sprintf(szBuffer,"FRAG 17  %.4f %.4f %.4f %.4f %.4f %.4f\n", FDB[i].a, FDB[i].b, FDB[i].c, FDB[i].al, FDB[i].be, FDB[i].ga);
        ar.WriteString(szBuffer);
        for(int j=0;j<FDB[i].atoms.size();j++) {
            Atom *tmpAt = &(FDB[i].atoms[j]);
            if(tmpAt->AtomType.CompareNoCase("Q")==0) {
                sprintf(szBuffer,"%s  %d   %.5f  %.5f  %.5f \n", tmpAt->AtomLabel, tmpAt->AtomicNumber, tmpAt->X, tmpAt->Y, tmpAt->Z);
            } else {
                sprintf(szBuffer,"%s  -%d   %.5f  %.5f  %.5f \n", tmpAt->AtomLabel, tmpAt->AtomicNumber, tmpAt->X, tmpAt->Y, tmpAt->Z);
            }
            ar.WriteString(szBuffer);
        }
        ar.WriteString("</"+FDB[i].tag_name+">\n\n"); 
    }    

	ar.Close();
	f.Close();
	return true;
}
bool IO::ReadFragmentDBFile(CString filename, vector<FragmentDBStruct> &FDB)
{
    CFileStatus status;
    if( CFile::GetStatus( filename, status ) == 0) {
        return false;
    }
  	CFile f(filename, CFile::modeRead);
	CArchive ar(&f,CArchive::load);
    CString line, line_upper;
    char szString[255];
    
    if(!ar.ReadString(line)) return false;
    
    while (ar.ReadString(line)) {
        if(line.Find("#")==0) {
            continue;
        } else if(line.Find("<")==0) {
            //new tag found
            FragmentDBStruct cfdb;
            line.Replace("<", "");
            line.Replace(">", "");
            cfdb.tag_name = line;
            while (ar.ReadString(line)) {
                line_upper = line;
                line_upper.MakeUpper();
                if(line_upper.Find("REM")==0) {
                    if(line_upper.Find("NAME:")!=-1) {
                        int q = line_upper.Find("NAME:");
                        cfdb.alt_names = line.Mid(q+6);
                    } else if(line_upper.Find("SRC:")!=-1) {
                        int q = line_upper.Find("SRC:");
                        cfdb.src = line.Mid(q+5);
                    } else {
                        cfdb.rem += line;
                        cfdb.rem += "\n";
                    }
                } else if(line_upper.Find("RESI")==0) {
                    cfdb.resi = line;
                } else if(line_upper.Find("FRAG")==0) {
                    int nb;
                    sscanf(line,"%s %d %f %f %f %f %f %f", &szString, &nb, &cfdb.a, &cfdb.b, &cfdb.c, &cfdb.al, &cfdb.be, &cfdb.ga);
                    while (ar.ReadString(line)) {
                        if(line.Find("</")==0) {
                            break;
                        } else {
                            float x, y, z;
                            int q = sscanf(line,"%s %d %f %f %f", &szString, &nb, &x, &y, &z);
                            if(q==5) {
                                /*
                                # e.g.  O1     1       x  y  z     <- Either the atom type is recognized by the atom name for 
                                #                                     positive numbers.
                                # e.g.  O1    -8       x  y  z     <- Or the atom type is defined by the negative atomic     
                                #                                     number.

                                */
                                if(nb<0) {
                                    nb = -1*nb;
                                    Atom tmpat(Atom::getAtomType(nb-1), szString, x, y, z);
                                    tmpat.AtomicNumber = nb;
                                    cfdb.atoms.push_back(tmpat);
                                } else {
                                    Atom tmpat(Atom::getAtomType(szString), szString, x, y, z);
                                    //tmpat.AtomicNumber = nb;
                                    cfdb.atoms.push_back(tmpat);
                                }
                            }
                        }

                    }
                } else if(line_upper.Find("JANA_RESTR")==0) {                    
                    int q = line_upper.Find("JANA_RESTR");
                    cfdb.jana_restr += line.Mid(q+11);
                    cfdb.jana_restr += "\n";
                } else if(line_upper.Find("CRYSTALS_RESTR")==0) {
                    int q = line_upper.Find("CRYSTALS_RESTR");
                    cfdb.Crystals_restr += line.Mid(q+15);
                    cfdb.Crystals_restr += "\n";
                } else {
                    if(line.Compare("")!=0) { //skip empty lines
                        cfdb.shelx_restr += line;
                        cfdb.shelx_restr += "\n";
                    }
                }
                if(line.Find("</")==0) {
                    cfdb.mol = getMoleculeFromFragmentDBEntry(cfdb);
                    FDB.push_back(cfdb);
                    break;
                }
            }
        }
    };

    ar.Close();
	f.Close();
    return true;
}
bool IO::ReadMolFile(CString filename, Molecule *mol)
{
/*
 1   benzene
 2   ACD/Labs0812062058
 3   
 4    6  6  0  0  0  0  0  0  0  0  1 V2000
 5      1.9050   -0.7932    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
 6      1.9050   -2.1232    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
 7      0.7531   -0.1282    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
 8      0.7531   -2.7882    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
 9     -0.3987   -0.7932    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
10     -0.3987   -2.1232    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
11    2  1  1  0  0  0  0
12    3  1  2  0  0  0  0
13    4  2  2  0  0  0  0
14    5  3  1  0  0  0  0
15    6  4  1  0  0  0  0
16    6  5  2  0  0  0  0
17   M  END
18   $$$$
*/   
    int nbAt;
    if(mol==0) return false;

    CString strInputData;
	CFile f(filename, CFile::modeRead);
	CArchive ar(&f,CArchive::load);

    //reading header
    if(!ar.ReadString(strInputData)) return false;
    if(!ar.ReadString(strInputData)) return false;
    if(!ar.ReadString(strInputData)) return false;
    
    if(!ar.ReadString(strInputData)) {
        return false;
    } else {
        //6  6  0  0  0  0  0  0  0  0  1 V2000
        vector<CString> values;
        if(!getValuesSeparatedByWhiteSpaces(strInputData, values)) return false;
        nbAt = atoi(values[0]);
    }
    //reading atoms
    for(int i=0;i<nbAt;i++) {
        if(!ar.ReadString(strInputData)) return false;
        vector<CString> values;
        if(!getValuesSeparatedByWhiteSpaces(strInputData, values)) return false;
        if(values.size()<4) return false;        
        mol->AddNewAtom(values[3], atof(values[0]), atof(values[1]), atof(values[2]));
    }	
    
    ar.Close();
	f.Close();

    return true;
}
bool IO::ReadCifFile(CString filename, UnitCell *pUnit)
{    
    if(pUnit==0) return false;

    CifReader cif;

    if(!cif.readFile(filename.GetString())) {
        AfxMessageBox("Error during file reading (" + filename + ")!"); 
        return false;
    }

    cif_crystal crystal = cif.getCrystal();
    std::vector<cif_atom> atoms = cif.getAtoms();

    pUnit->A = crystal.a;
    pUnit->B = crystal.b;
    pUnit->C = crystal.c;
    pUnit->Alpha = crystal.al;
    pUnit->Beta = crystal.be;
    pUnit->Gama = crystal.ga;
    
    //std::cout<<"matrixes"<<std::endl;
    if(crystal.equiv_pos.size()!=0) {
        pUnit->SetSpaceGroup(crystal.sg.c_str(), false);
        pUnit->listMatrix = crystal.equiv_pos_matrixes;        
    } else {
        pUnit->SetSpaceGroup(crystal.sg.c_str(), true);
    }

    for(int i=0;i<atoms.size();i++) {
        pUnit->AddNewAtomInFract(atoms[i].type.c_str(), atoms[i].name.c_str(), atoms[i].x, atoms[i].y, atoms[i].z);
    }

    return true;
}
bool IO::ReadAtomFile(CString filename, UnitCell *pUnit)
{
	float x, y, z, r, t, occ, xxx;
	int nb=-1;
    if(pUnit==0) return false;

    CString strInputData;
    char szString[100];
	CFile f(filename, CFile::modeRead);
	CArchive ar(&f,CArchive::load);
	
    int nbReads = 0;
	//nbAtoms   nbMolecules   rotMatrix
	//3          0              0                0
    //reading nb of atom of all phases
    //saving nb of atoms of the first phase only
    do { 
	    if(!ar.ReadString(strInputData)) return false;
        int tnb=0;
        nbReads = sscanf(strInputData,"%d %f %f %f %f", &tnb, &x, &y, &z, &xxx);
        if(nb==-1) nb = tnb;
    } while(nbReads==4);
     
	//if(!ar.ReadString(strInputData)) return false;
	if(!ar.ReadString(strInputData)) return false;
	if(!ar.ReadString(strInputData)) return false;
	if(!ar.ReadString(strInputData)) return false;
	
	//Cu1       1  1     0.000000 0.347375 0.250000 0.258493
	for(int i=0;i<nb;i++) {
		if(!ar.ReadString(strInputData)) return false;
        sscanf(strInputData,"%s %d %d %f %f %f %f", szString, &t, &r, &occ, &x, &y, &z);        
		pUnit->AddNewAtomInFract(Atom::getAtomType(szString), szString, x, y, z);
		if(!ar.ReadString(strInputData)) return false;
	}
    
    ar.Close();
	f.Close();

    return true;
}
bool IO::ReadCavityFile(CString filename, UnitCell *pUnit)
{
    float x, y, z, r, t;
    if(pUnit==0) return false;

    CString strInputData;
    char szString[100];
	CFile f(filename, CFile::modeRead);
	CArchive ar(&f,CArchive::load);

	do{
		if(!ar.ReadString(strInputData)) break;
        sscanf(strInputData,"%f	%f %f %f %f %s", &x, &y, &z, &r, &t, szString);
        pUnit->AddNewCavity(x, y, z, r, t, szString);
	}while(true);
    
    ar.Close();
	f.Close();

    return true;
}
void IO::GetFirstParamFromString(CString &s, CString &lparam, CString &rparam)
{
   s.TrimLeft();
   int i = s.Find(' ', 0);
   if(i==-1) {
       lparam=s;
       rparam="";
       return;
   }
   lparam = s.Mid(0, i);
   rparam = s.Mid(i+1);     
}
bool IO::ReadMCEFile(CString filename, UnitCell *pUnit[], VoxelMap *pMap[], 
                     DisplayStyleStruct &DSS, SpaceGroupSetup &SGS, LevelControlSetup &LCS, vector<Atom> &PL, SelectStruct &LS,
					 float *curquat, float *lastquat, float *shift, float &scale)
{
    //if(pUnit[0]==NULL) return false;
    //if(pMap[0]==NULL) return false;    
    CString strInputData;
    map<CString, CString> parameters, *currmappar;
    map<CString, CString> mappar[2];
    CFile f;
	if(!f.Open(filename, CFile::modeRead)) {
		AfxMessageBox("File: " + filename + " does not exists."); 
		return false;
	}
	CArchive ar(&f,CArchive::load);
    CString param, value;
    VoxelMap *currMap = NULL;
    vector<Atom> atoms;
    do{
		if(!ar.ReadString(strInputData)) break;
        GetFirstParamFromString(strInputData, param, value);
        if(param.Compare("listMatrix")==0) {
            if(pUnit[0]==NULL) pUnit[0] = new UnitCell();
            //save operations of symmetry
            vector<vector<float>> tmp;
            tmp.resize(3);
            tmp[0].resize(4);
            tmp[1].resize(4);
            tmp[2].resize(4);
            sscanf(value, "%f %f %f %f %f %f %f %f %f %f %f %f", &tmp[0][0], &tmp[0][1], &tmp[0][2], &tmp[0][3],
                                                                 &tmp[1][0], &tmp[1][1], &tmp[1][2], &tmp[1][3],
                                                                 &tmp[2][0], &tmp[2][1], &tmp[2][2], &tmp[2][3]);   
            pUnit[0]->listMatrix.push_back(tmp);
        } else if(param.Compare("m_original_atoms_fract")==0) {
            if(pUnit[0]==NULL) pUnit[0] = new UnitCell();
            char label[100], type[100];
            float x, y, z;
            int k;
            sscanf(value, "%s %s %d %f %f %f", label, type, &k, &x, &y, &z);           
            atoms.push_back(Atom(type, label, x, y, z));
            //pUnit[0]->AddNewAtomInFract(type, label, x, y, z);
        } else if(param.Compare("peak_atom")==0) {
            char label[100], type[100];
            float x, y, z;
            int k;
            sscanf(value, "%s %s %f %f %f", label, type, &x, &y, &z);
            PL.push_back(Atom(type, label, x, y, z));
        } else if(param.Compare("Map_object_id")==0) {
            parameters[param] = value;
            //define the current map
            int nb;
            sscanf(value, "%d", &nb);
            if((nb>1) || (nb<0)) {
                AfxMessageBox("Error: Map_object_id = " + value + " is has not been expected!\n Set Map_object_id from 0 to 1!");
                return false;
            }
            if(pMap[nb]==NULL) pMap[nb] = new VoxelMap();
            currMap = pMap[nb];
            currmappar = &mappar[nb];
        } else if(param.Compare("Voxels:")==0) {
            if(currMap==NULL) {
                AfxMessageBox("Error: Voxels found, but Map_object_id is not defined!\n Define Map_object_id before Voxels!");
                return false;
            }
            float val;
            sscanf(value, "%f", &val);
            if(((*currmappar).count("Lines")==0) || ((*currmappar).count("Columns")==0) || ((*currmappar).count("Sections")==0)) {
                AfxMessageBox("Error: Define Lines, Columns and Sections before you define Voxels!");
                return false;
            }
            currMap->Lines = stringToInt((*currmappar)["Lines"].GetString());
            currMap->Columns = stringToInt((*currmappar)["Columns"].GetString());
            currMap->Sections = stringToInt((*currmappar)["Sections"].GetString());
            
            currMap->Voxel.clear();
            currMap->MemoryForVoxel();                       
            for(int i=0;i<currMap->Voxel.size();i++) {                
                for(int j=0;j<currMap->Voxel[i].size();j++) {
                    for(int k=0;k<currMap->Voxel[i][j].size();k++) {
                        if(!ar.ReadString(strInputData)) {
                            AfxMessageBox("Error: Unexpected end of file while reading Voxels!");
                            return false;
                        }
                        currMap->Voxel[i][j][k] = atof(strInputData);
                    }
                }
            }
        } else {
            if(parameters.count("Map_object_id")>0) {
                (*currmappar)[param]=value;
            } else {
                parameters[param] = value;
            }
        }
	} while(true);

    LoadDisplayStyle(parameters, &DSS);
    LoadSGSetup(parameters, &SGS);
    LoadLabelStruct(parameters, &LS);
    
    if(!parameters["curquat"].IsEmpty()) {
        sscanf(parameters["curquat"], "%f %f %f %f", &curquat[0], &curquat[1], &curquat[2], &curquat[3]);                                                                                      
    }
    if(!parameters["lastquat"].IsEmpty()) {
        sscanf(parameters["lastquat"], "%f %f %f %f", &lastquat[0], &lastquat[1], &lastquat[2], &lastquat[3]);                                                                                      
    }
    if(!parameters["shift"].IsEmpty()) {
        sscanf(parameters["shift"], "%f %f %f", &shift[0], &shift[1], &shift[2]);                                                                                      
    }
    if(!parameters["scale"].IsEmpty()) {
        scale = atof(parameters["scale"]);
    }
    if(!parameters["iMapMode"].IsEmpty()) {
        sscanf(parameters["iMapMode"], "%d %d", &LCS.iMapMode[0],&LCS.iMapMode[1]);
    }
    if(!parameters["iScrollBarMode"].IsEmpty()) {
        sscanf(parameters["iScrollBarMode"], "%d %d", &LCS.iScrollBarMode[0],&LCS.iScrollBarMode[1]);
    }
    /*  

    //levelcontrols
    sprintf(szBuffer,"iMapMode %d %d\n", LCS.iMapMode[0], LCS.iMapMode[1]);
    ar.WriteString(szBuffer);
    sprintf(szBuffer,"iScrollBarMode %d %d\n", LCS.iScrollBarMode[0], LCS.iScrollBarMode[1]);
    ar.WriteString(szBuffer);
    */

    if(pMap[0]!=NULL) {
        pMap[0]->LoadValues(mappar[0]);
    }
    if(pMap[1]!=NULL) {
        pMap[1]->LoadValues(mappar[1]);
    }
    if(pUnit[0]!=NULL) {
        pUnit[0]->LoadValues(parameters);
        for(int i=0;i<atoms.size();i++) {
            pUnit[0]->AddNewAtomInFract(atoms[i]);
        }
        pUnit[0]->generateMoleculesFromAssymPart();
    }

    ar.Close();
	f.Close();
    return true;
}
 bool IO::SaveMCEInput(     CString filename, UnitCell *pUnit[], VoxelMap *pMap[], 
							DisplayStyleStruct DSS, 
							vector<ElementDescription> ED, 
							SpaceGroupSetup SGS,
							LevelControlSetup LCS,
                            vector<Atom> PL,
                            SelectStruct LS,
							float curquat[4], float lastquat[4], float shift[3], float scale)
 {        
	if(pUnit[0]==0) return false;
	if(pMap[0]==0) return false;
    char szBuffer[255];

    //LZOpenFile


    CFile f( filename, CFile::modeCreate | CFile::modeWrite);
    CArchive ar(&f,CArchive::store);
    
    //save view
    sprintf(szBuffer,"curquat %f %f %f %f\n", curquat[0], curquat[1], curquat[2], curquat[3]);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"lastquat %f %f %f %f\n", lastquat[0], lastquat[1], lastquat[2], lastquat[3]);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"shift %f %f %f\n", shift[0], shift[1], shift[2]);
	ar.WriteString(szBuffer);
    sprintf(szBuffer,"scale %f\n", scale);
    ar.WriteString(szBuffer);
    
    //save display style
    ar.WriteString(SaveDisplayStyle(DSS));
    SaveLabelStruct(LS, &ar);
    
    //save differencies in atomic types
    for(int i=0;i<ED.size();i++) {
        sprintf(szBuffer,"Atomic_type %s %f %f %f %f %f %f %f\n", ED[i].Name, ED[i].ATWE, ED[i].COVR, ED[i].T, ED[i].VDWR, ED[i].color[0], ED[i].color[1], ED[i].color[2]);
    }



    //save spacegroup and unit cell limits
    sprintf(szBuffer,"m_bSpaceGroup %s\n", SGS.m_bSpaceGroup);
    ar.WriteString(szBuffer);
    sprintf(szBuffer,"m_bvisible %d\n", SGS.m_bvisible);
    ar.WriteString(szBuffer);
    sprintf(szBuffer,"m_bFromX_Axis %d\n", SGS.m_bFromX_Axis);
    ar.WriteString(szBuffer);
    sprintf(szBuffer,"m_bFromY_Axis %d\n", SGS.m_bFromY_Axis);
    ar.WriteString(szBuffer);
    sprintf(szBuffer,"m_bFromZ_Axis %d\n", SGS.m_bFromZ_Axis);
    ar.WriteString(szBuffer);
    sprintf(szBuffer,"m_bToX_Axis %d\n", SGS.m_bToX_Axis);
    ar.WriteString(szBuffer);
    sprintf(szBuffer,"m_bToY_Axis %d\n", SGS.m_bToY_Axis);
    ar.WriteString(szBuffer);
    sprintf(szBuffer,"m_bToZ_Axis %d\n", SGS.m_bToZ_Axis);
    ar.WriteString(szBuffer);

    //levelcontrols
    sprintf(szBuffer,"iMapMode %d %d\n", LCS.iMapMode[0], LCS.iMapMode[1]);
    ar.WriteString(szBuffer);
    sprintf(szBuffer,"iScrollBarMode %d %d\n", LCS.iScrollBarMode[0], LCS.iScrollBarMode[1]);
    ar.WriteString(szBuffer);
    
	for(int i=0;i<PL.size();i++) {
        sprintf(szBuffer,"peak_atom %s %s %f %f %f\n", PL[i].AtomLabel, PL[i].AtomType, PL[i].X, PL[i].Y, PL[i].Z);
        ar.WriteString(szBuffer);
	}
        
    //Unit Cell
    pUnit[0]->SaveToFile(&ar);

    for(int i=0;i<NMAPS;i++) {
        if(pMap[i]==NULL) continue;
        sprintf(szBuffer,"Map_object_id %d\n", i);
        ar.WriteString(szBuffer);  
        pMap[i]->SaveToFile(&ar);
    }

    ar.Close();
	f.Close();
    return true;
 }
 void IO::LoadSGSetup(map<CString, CString> params, SpaceGroupSetup *res)
 {
     if(!params["m_bSpaceGroup"].IsEmpty()) {
        res->m_bSpaceGroup = params["m_bSpaceGroup"];
     }
     if(!params["m_bvisible"].IsEmpty()) {
        res->m_bvisible = atoi(params["m_bvisible"]);
     }
     if(!params["m_bFromX_Axis"].IsEmpty()) {
        res->m_bFromX_Axis = atof(params["m_bFromX_Axis"]);
     }
     if(!params["m_bFromY_Axis"].IsEmpty()) {
        res->m_bFromY_Axis = atof(params["m_bFromY_Axis"]);
     }
     if(!params["m_bFromZ_Axis"].IsEmpty()) {
        res->m_bFromZ_Axis = atof(params["m_bFromZ_Axis"]);
     }
     if(!params["m_bToX_Axis"].IsEmpty()) {
        res->m_bToX_Axis = atof(params["m_bToX_Axis"]);
     }
     if(!params["m_bToY_Axis"].IsEmpty()) {
        res->m_bToY_Axis = atof(params["m_bToY_Axis"]);
     }
     if(!params["m_bToZ_Axis"].IsEmpty()) {
        res->m_bToZ_Axis = atof(params["m_bToZ_Axis"]);
     }
 }
 void IO::LoadLabelStruct(map<CString, CString> params, SelectStruct *res)
 {
     if(!params["m_bMolecule"].IsEmpty()) {
        res->m_bMolecule = atoi(params["m_bMolecule"]);
     }
     //if(!params["m_bNatoms"].IsEmpty()) {
     //   res->m_bNatoms = atoi(params["m_bNatoms"]);
     //}
     //if(!params["m_bBonding"].IsEmpty()) {
     //   res->m_bBonding = atoi(params["m_bBonding"]);
     //}
     //if(!params["m_bHits"].IsEmpty()) {
     //   res->m_bHits = atoi(params["m_bHits"]);
     //}
     if(!params["m_rotate_selected_molecule"].IsEmpty()) {
        res->m_rotate_selected_molecule = atoi(params["m_rotate_selected_molecule"]);
     }
     if(!params["m_use_atom_instead_of_mass_center"].IsEmpty()) {
        res->m_use_atom_instead_of_mass_center = atoi(params["m_use_atom_instead_of_mass_center"]);
     }
     if(!params["m_track_ball_rot_m"].IsEmpty()) {
        sscanf(params["m_track_ball_rot_m"], "%f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f", &res->m_track_ball_rot_m[0][0], &res->m_track_ball_rot_m[0][1], &res->m_track_ball_rot_m[0][2], &res->m_track_ball_rot_m[0][3],
                                                                                                &res->m_track_ball_rot_m[1][0], &res->m_track_ball_rot_m[1][1], &res->m_track_ball_rot_m[1][2], &res->m_track_ball_rot_m[1][3],
                                                                                                &res->m_track_ball_rot_m[2][0], &res->m_track_ball_rot_m[2][1], &res->m_track_ball_rot_m[2][2], &res->m_track_ball_rot_m[2][3],
                                                                                                &res->m_track_ball_rot_m[3][0], &res->m_track_ball_rot_m[3][1], &res->m_track_ball_rot_m[3][2], &res->m_track_ball_rot_m[3][3]);
     }
     if(!params["m_bName"].IsEmpty()) {
        res->m_bName = params["m_bName"];
     }
     if(!params["position"].IsEmpty()) {
        sscanf(params["position"], "%f %f %f ", &res->position[0], &res->position[1], &res->position[2]);
     }
     if(!params["positionID"].IsEmpty()) {
        res->positionID = atoi(params["positionID"]);
     }
     if(!params["Rot"].IsEmpty()) {
        sscanf(params["Rot"], "%f %f %f %f %f %f %f %f %f", &res->Rot[0][0], &res->Rot[0][1], &res->Rot[0][2],
                                                            &res->Rot[1][0], &res->Rot[1][1], &res->Rot[1][2],
                                                            &res->Rot[2][0], &res->Rot[2][1], &res->Rot[2][2]);
     }
     if(!params["Trans"].IsEmpty()) {
        sscanf(params["Trans"], "%f %f %f", &res->Trans[0], &res->Trans[1], &res->Trans[2]);
     }
     if(!params["m_apply_symmetry"].IsEmpty()) {
        res->m_apply_symmetry = atoi(params["m_apply_symmetry"]);
     }     
     if(!params["m_move_fragment_to_mouse_click"].IsEmpty()) {
        res->m_move_fragment_to_mouse_click = atoi(params["m_move_fragment_to_mouse_click"]);
     }
     if(!params["m_is_selected_atom_part_of_selected_fragment"].IsEmpty()) {
        res->m_is_selected_atom_part_of_selected_fragment = atoi(params["m_is_selected_atom_part_of_selected_fragment"]);
     }     

 }
 void IO::SaveLabelStruct(SelectStruct LS, CArchive *ar)
 {
     char szBuffer[255];
     /*
	int m_bMolecule;
	int m_bNatoms;
	int m_bBonding;
	int m_bHits;
    bool m_rotate_selected_molecule;
    bool m_use_atom_instead_of_mass_center;
    float m_track_ball_rot_m[4][4];
	CString m_bName;
    float position[3];
    Atom *pAtom;
    Atom *pPeak;
    int positionID;
	float Rot[3][3];
	float Trans[3];
    bool m_apply_symmetry;
    bool m_move_fragment_to_mouse_click;
    bool m_is_selected_atom_part_of_selected_fragment;

     */
    sprintf(szBuffer,"m_bMolecule %d\n", LS.m_bMolecule);
    ar->WriteString(szBuffer);
    //sprintf(szBuffer,"m_bNatoms %d\n", LS.m_bNatoms);
    //ar->WriteString(szBuffer);
    //sprintf(szBuffer,"m_bBonding %d\n", LS.m_bBonding);
    //ar->WriteString(szBuffer);
    //sprintf(szBuffer,"m_bHits %d\n", LS.m_bHits);
    //ar->WriteString(szBuffer);
    sprintf(szBuffer,"m_rotate_selected_molecule %d\n", LS.m_rotate_selected_molecule);
    ar->WriteString(szBuffer);
    sprintf(szBuffer,"m_use_atom_instead_of_mass_center %d\n", LS.m_use_atom_instead_of_mass_center);
    ar->WriteString(szBuffer);
    sprintf(szBuffer,"m_track_ball_rot_m %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n", LS.m_track_ball_rot_m[0][0], LS.m_track_ball_rot_m[0][1], LS.m_track_ball_rot_m[0][2], LS.m_track_ball_rot_m[0][3],
                                                                                             LS.m_track_ball_rot_m[1][0], LS.m_track_ball_rot_m[1][1], LS.m_track_ball_rot_m[1][2], LS.m_track_ball_rot_m[1][3],
                                                                                             LS.m_track_ball_rot_m[2][0], LS.m_track_ball_rot_m[2][1], LS.m_track_ball_rot_m[2][2], LS.m_track_ball_rot_m[2][3],
                                                                                             LS.m_track_ball_rot_m[3][0], LS.m_track_ball_rot_m[3][1], LS.m_track_ball_rot_m[3][2], LS.m_track_ball_rot_m[3][3]);
    ar->WriteString(szBuffer);
    sprintf(szBuffer,"m_bName %s\n", LS.m_bName);
    ar->WriteString(szBuffer);
    sprintf(szBuffer,"position %f %f %f\n", LS.position[0], LS.position[1], LS.position[2]);
    ar->WriteString(szBuffer);
    sprintf(szBuffer,"positionID %d\n", LS.positionID);
    ar->WriteString(szBuffer);
    sprintf(szBuffer,"Rot %f %f %f %f %f %f %f %f %f\n", LS.Rot[0][0],LS.Rot[0][1],LS.Rot[0][2],
                                                         LS.Rot[1][0],LS.Rot[1][1],LS.Rot[1][2],
                                                         LS.Rot[2][0],LS.Rot[2][1],LS.Rot[2][2]);
    ar->WriteString(szBuffer);
    sprintf(szBuffer,"Trans %f %f %f\n", LS.Trans[0], LS.Trans[1], LS.Trans[2]);
    ar->WriteString(szBuffer);
    sprintf(szBuffer,"m_apply_symmetry %d\n", LS.m_apply_symmetry);
    ar->WriteString(szBuffer);
    sprintf(szBuffer,"m_move_fragment_to_mouse_click %d\n", LS.m_move_fragment_to_mouse_click);
    ar->WriteString(szBuffer);
    sprintf(szBuffer,"m_is_selected_atom_part_of_selected_fragment %d\n", LS.m_is_selected_atom_part_of_selected_fragment);
    ar->WriteString(szBuffer);

 }
  void IO::LoadDisplayStyle(map<CString, CString> params, DisplayStyleStruct *res)
 {

    if(!params["DisplayStyle"].IsEmpty()) {
        res->DisplayStyle = (DisplayStyleMode) atoi(params["DisplayStyle"]);
    }
    if(!params["MapDisplayStyle"].IsEmpty()) {
        res->MapDisplayStyle = (MapDisplayStyleMode) atoi(params["MapDisplayStyle"]);
    }
    if(!params["iBondsColorStyle"].IsEmpty()) {
        res->iBondsColorStyle = (BondsColorStyle) atoi(params["iBondsColorStyle"]);
    }
    if(!params["LineSize"].IsEmpty()) {
        res->LineSize = atof(params["LineSize"]);
    }
    if(!params["BallSize"].IsEmpty()) {
        res->BallSize = atof(params["BallSize"]);
    }
    if(!params["ScaledFactor"].IsEmpty()) {
        res->ScaledFactor = atof(params["ScaledFactor"]);
    }
    if(!params["fCylinderSize"].IsEmpty()) {
        res->fCylinderSize = atof(params["fCylinderSize"]);
    }
    if(!params["fScaledBallResolution"].IsEmpty()) {
        res->fScaledBallResolution = atof(params["fScaledBallResolution"]);
    }
    if(!params["fCylinderResolution"].IsEmpty()) {
        res->fCylinderResolution = atof(params["fCylinderResolution"]);
    }
    if(!params["fBallResolution"].IsEmpty()) {
        res->fBallResolution = atof(params["fBallResolution"]);
    }
    if(!params["fMapLineSize"].IsEmpty()) {
        res->fMapLineSize = atof(params["fMapLineSize"]);
    }
    if(!params["fMapCylinderSize"].IsEmpty()) {
        res->fMapCylinderSize = atof(params["fMapCylinderSize"]);
    }
    if(!params["fUnitCellLineSize"].IsEmpty()) {
        res->fUnitCellLineSize = atof(params["fUnitCellLineSize"]);
    }
    if(!params["fMapCylinderResolution"].IsEmpty()) {
        res->fMapCylinderResolution = atof(params["fMapCylinderResolution"]);
    }
    if(!params["BondsColor"].IsEmpty()) {
        sscanf(params["BondsColor"], "%f %f %f", &res->BondsColor[0], &res->BondsColor[1], &res->BondsColor[2]);
    }
    if(!params["LabelsColor"].IsEmpty()) {
        sscanf(params["LabelsColor"], "%f %f %f", &res->LabelsColor[0], &res->LabelsColor[1], &res->LabelsColor[2]);
    }
    if(!params["BgColor"].IsEmpty()) {
        sscanf(params["BgColor"], "%f %f %f", &res->BgColor[0], &res->BgColor[1], &res->BgColor[2]);
    }
    if(!params["iLabelSize"].IsEmpty()) {
        res->iLabelSize = atoi(params["iLabelSize"]);
    }
    if(!params["bDrawMapEdges"].IsEmpty()) {
        res->bDrawMapEdges = atoi(params["bDrawMapEdges"]);
    }
    if(!params["bUnitCellVisible"].IsEmpty()) {
        res->bUnitCellVisible = atoi(params["bUnitCellVisible"]);
    }
    if(!params["bBondsVisible"].IsEmpty()) {
        res->bBondsVisible = atoi(params["bBondsVisible"]);
    }
    if(!params["bOrganic"].IsEmpty()) {
        res->bOrganic = atoi(params["bOrganic"]);
    }
    if(!params["bLoadMapSettings"].IsEmpty()) {
        res->bLoadMapSettings = atoi(params["bLoadMapSettings"]);
    }
    if(!params["bShowLabels"].IsEmpty()) {
        res->bShowLabels = atoi(params["bShowLabels"]);
    }
    if(!params["bShowHLabels"].IsEmpty()) {
        res->bShowHLabels = atoi(params["bShowHLabels"]);
    }
    if(!params["bShowPeaksCloseContacts"].IsEmpty()) {
        res->bShowPeaksCloseContacts = atoi(params["bShowPeaksCloseContacts"]);
    }
    if(!params["bShowPeaksSymmetryEquivalents"].IsEmpty()) {
        res->bShowPeaksSymmetryEquivalents = atoi(params["bShowPeaksSymmetryEquivalents"]);
    }
 }
 CString IO::SaveDisplayStyle(DisplayStyleStruct DSS)
 {
     CString res;
     char szBuffer[255];
     /*
     DisplayStyleMode DisplayStyle;
	MapDisplayStyleMode	MapDisplayStyle;
	BondsColorStyle iBondsColorStyle;
	float LineSize;
	float BallSize;
	float ScaledFactor;
	float fCylinderSize;
	float fScaledBallResolution;
	float fCylinderResolution;
	float fBallResolution;
	float fMapLineSize;
	float fMapCylinderSize;
	float fUnitCellLineSize;
	float fMapCylinderResolution;
	float BondsColor[3];//R,G,B
    float LabelsColor[3];    
	float MapsColors[NMAPS][3][3];
    float MapsMSHColors[NMAPS][2][3];
	float BgColor[3];
    int  iLabelSize;
    BOOL bDrawMapEdges;
	BOOL bUnitCellVisible;
	BOOL bBondsVisible;
	BOOL bOrganic;
	BOOL bLoadMapSettings;
    BOOL bShowLabels;
    BOOL bShowHLabels;
    BOOL bShowPeaksCloseContacts;
    BOOL bShowPeaksSymmetryEquivalents;
     */
    
    sprintf(szBuffer,"DisplayStyle %d\n", DSS.DisplayStyle);
    res+=szBuffer;
    sprintf(szBuffer,"MapDisplayStyle %d\n", DSS.MapDisplayStyle);
    res+=szBuffer;
    sprintf(szBuffer,"iBondsColorStyle %d\n", DSS.iBondsColorStyle);
    res+=szBuffer;
    sprintf(szBuffer,"LineSize %f\n", DSS.LineSize);
    res+=szBuffer;
    sprintf(szBuffer,"BallSize %f\n", DSS.BallSize);
    res+=szBuffer;
    sprintf(szBuffer,"ScaledFactor %f\n", DSS.ScaledFactor);
    res+=szBuffer;
    sprintf(szBuffer,"fCylinderSize %f\n", DSS.fCylinderSize);
    res+=szBuffer;
    sprintf(szBuffer,"fScaledBallResolution %f\n", DSS.fScaledBallResolution);
    res+=szBuffer;
    sprintf(szBuffer,"fCylinderResolution %f\n", DSS.fCylinderResolution);
    res+=szBuffer;
    sprintf(szBuffer,"fBallResolution %f\n", DSS.fBallResolution);
    res+=szBuffer;
    sprintf(szBuffer,"fMapLineSize %f\n", DSS.fMapLineSize);
    res+=szBuffer;
    sprintf(szBuffer,"fMapCylinderSize %f\n", DSS.fMapCylinderSize);
    res+=szBuffer;
    sprintf(szBuffer,"fUnitCellLineSize %f\n", DSS.fUnitCellLineSize);
    res+=szBuffer;
    sprintf(szBuffer,"fMapCylinderResolution %f\n", DSS.fMapCylinderResolution);
    res+=szBuffer;
    sprintf(szBuffer,"BondsColor %f %f %f\n", DSS.BondsColor[0], DSS.BondsColor[1], DSS.BondsColor[2]);
    res+=szBuffer;
    /*
    sprintf(szBuffer,"MapsColors %f %f %f %f %f %f %f %f %f\n", DSS.MapsColors[0][0], DSS.MapsColors[0][1], DSS.MapsColors[0][2],
                                                                  DSS.MapsColors[1][0], DSS.MapsColors[1][1], DSS.MapsColors[1][2], 
                                                                  DSS.MapsColors[2][0], DSS.MapsColors[2][1], DSS.MapsColors[2][2]);
    */
    /*
    res+=szBuffer;
    sprintf(szBuffer,"MapsMSHColors %f %f %f %f %f %f\n", DSS.MapsMSHColors[0][0], DSS.MapsMSHColors[0][1], DSS.MapsMSHColors[0][2], 
                                                            DSS.MapsMSHColors[1][0], DSS.MapsMSHColors[1][1], DSS.MapsMSHColors[1][2]);
    */
    res+=szBuffer;
    sprintf(szBuffer,"BgColor %f %f %f\n", DSS.BgColor[0], DSS.BgColor[1], DSS.BgColor[2]);
    res+=szBuffer;
    sprintf(szBuffer,"iLabelSize %d\n", DSS.iLabelSize);
    res+=szBuffer;
    sprintf(szBuffer,"bDrawMapEdges %d\n", DSS.bDrawMapEdges);
    res+=szBuffer;
    sprintf(szBuffer,"bUnitCellVisible %d\n", DSS.bUnitCellVisible);
    res+=szBuffer;
    sprintf(szBuffer,"bBondsVisible %d\n", DSS.bBondsVisible);
    res+=szBuffer;
    sprintf(szBuffer,"bOrganic %d\n", DSS.bOrganic);
    res+=szBuffer;
    sprintf(szBuffer,"bLoadMapSettings %d\n", DSS.bLoadMapSettings);
    res+=szBuffer;
    sprintf(szBuffer,"bShowLabels %d\n", DSS.bShowLabels);
    res+=szBuffer;
    sprintf(szBuffer,"bShowHLabels %d\n", DSS.bShowHLabels);
    res+=szBuffer;
    sprintf(szBuffer,"bShowPeaksCloseContacts %d\n", DSS.bShowPeaksCloseContacts);
    res+=szBuffer;
    sprintf(szBuffer,"bShowPeaksSymmetryEquivalents %d\n", DSS.bShowPeaksSymmetryEquivalents);
    res+=szBuffer;

    return res;
}
std::string IO::floatToString(float n)
{
    std::ostringstream os("");
    os << n;
    return os.str();
}
std::string IO::floatToString(float n, int nb_decimal)
{
    std::ostringstream os("");
    os << std::fixed;
    os << std::setprecision(nb_decimal) << n; 
    return os.str();
}
float IO::stringToFloat(std::string str)
{
    std::istringstream is(str);
    double x;
    if (!(is >> x))
        return 0;
    return x;
}
float IO::stringToInt(std::string str)
{
    std::istringstream is(str);
    int x;
    if (!(is >> x))
        return 0;
    return x;
}
bool IO::SavePovRayInput(	CString filename, UnitCell *pUnit[], VoxelMap *pMap[], 
							DisplayStyleStruct DSS, 
							ElementDescription ED[255], 
							SpaceGroupSetup SGS,
							LevelControlSetup LCS,
							float matrix[16])
{
	int nbUnit = 0, nbMap = 0;
	
	for(int i=0;i<2;i++){
		if(pUnit[i]!=0) nbUnit++;
		if(pMap[i]!=0) nbMap++;
	}
	if((nbUnit==0)&&(nbMap==0)) return false;
	if((pUnit[0]==0)&&(nbUnit!=0)) return false;
	if((pMap[0]==0)&&(nbMap!=0)) return false;


	char szBuffer[250];
	CString szAtom;
	CFile f( filename, CFile::modeCreate | CFile::modeWrite);
	CArchive ar(&f,CArchive::store);

	// make a copy of the POV header file to the destination ...
	VERIFY(::GetModuleFileName(AfxGetInstanceHandle(), szBuffer, _MAX_PATH));

	CString sPath = (CString)szBuffer;
	sPath = sPath.Left(sPath.ReverseFind('\\'));	
	sPath=sPath+"\\MCE_POVray_head.def";
	
	CFile f_header(sPath, CFile::modeRead);
	CArchive ar_header(&f_header,CArchive::load);
	CString strInputData;

	ar_header.ReadString(strInputData);
	do
	{
		if(strInputData.Left(13)=="//END OF FILE") break;
		ar.WriteString(strInputData + "\n");
	}while(ar_header.ReadString(strInputData));

	ar_header.Close();
	f_header.Close();

	//background color
	sprintf(szBuffer,"background { color rgb <%f, %f,%f>}\n", DSS.BgColor[0], DSS.BgColor[1], DSS.BgColor[2]);
	ar.WriteString(szBuffer);

	sprintf(szBuffer,"#declare cellcylRadius=%0.3f;\n", DSS.fUnitCellLineSize/100);
	ar.WriteString(szBuffer);

	//define bonds
	if((DSS.DisplayStyle==Cylinder)|(DSS.DisplayStyle==Ball_and_cylinder)|(DSS.DisplayStyle==Scaled_ball_and_cylinder))
	    sprintf(szBuffer,"#declare bondcylRadius=%0.3f;", DSS.fCylinderSize);
	else
		sprintf(szBuffer,"#declare bondcylRadius=0.04;");//bonds definition -> Lines
		
	ar.WriteString(szBuffer);
	ar.WriteString("\n");

	if(DSS.iBondsColorStyle==OneColor) {
		sprintf(szBuffer,"#declare bondsColor = color rgb <%f, %f, %f>;\n", DSS.BondsColor[0], DSS.BondsColor[1], DSS.BondsColor[2]);//bonds definition -> Lines
		ar.WriteString(szBuffer);
	}
	
	//define Atoms
	if (DSS.DisplayStyle != Scaled_ball_and_cylinder)
	{
		if (DSS.DisplayStyle==Cylinder)
			sprintf(szBuffer,"#declare AtomRadius=%0.3f;", DSS.fCylinderSize);
		else
			sprintf(szBuffer,"#declare AtomRadius=%0.3f;", DSS.BallSize);
	}
	else
		sprintf(szBuffer,"#declare ScaledFactor=%f;", DSS.ScaledFactor);
	
	ar.WriteString(szBuffer);
	ar.WriteString("\n");
	
	//Map
	if (DSS.MapDisplayStyle == AsCylinder)
		sprintf(szBuffer,"#declare MapcylRadius=%f;", DSS.fMapCylinderSize);
	  else
		sprintf(szBuffer,"#declare MapcylRadius=0.015;");
	ar.WriteString(szBuffer);
	ar.WriteString("\n");
	//map colors
	for(int i=0;i<nbMap;i++){
		sprintf(szBuffer,"#declare Map%dLevel0Color= color rgb <%f, %f,%f>;\n", i, pMap[i]->layer_color[0][0], pMap[i]->layer_color[0][1], pMap[i]->layer_color[0][2]);
		ar.WriteString(szBuffer);
		sprintf(szBuffer,"#declare Map%dLevel1Color= color rgb <%f, %f,%f>;\n", i, pMap[i]->layer_color[1][0], pMap[i]->layer_color[1][1], pMap[i]->layer_color[1][2]);
		ar.WriteString(szBuffer);
		sprintf(szBuffer,"#declare Map%dLevel2Color= color rgb <%f, %f,%f>;\n", i, pMap[i]->layer_color[2][0], pMap[i]->layer_color[2][1], pMap[i]->layer_color[2][2]);
		ar.WriteString(szBuffer);
	}

	//element description
	for(int i=0;i<105;i++)
	{
		szAtom = ED[i].Name.MakeUpper();		
		sprintf(szBuffer,"#declare ColorFor%0.2s  = color rgb <%f, %f,%f>;\n", szAtom, ED[i].color[0], ED[i].color[1], ED[i].color[2]);				
		ar.WriteString(szBuffer);
		sprintf(szBuffer,"#declare AtomRadiusFor%0.2s  = %f;\n", szAtom, ED[i].VDWR);
		ar.WriteString(szBuffer);
	}

	// drawing of the scene
	ar.WriteString("union {\n");
	ar.WriteString("union {\n");
	
	//TODO: POVDrawMark(); 
	//POVDrawPeaks(&ar);

	for(int i=SGS.m_bFromX_Axis;i<SGS.m_bToX_Axis+1;i++)
		for(int j=SGS.m_bFromY_Axis;j<SGS.m_bToY_Axis+1;j++)
			for(int k=SGS.m_bFromZ_Axis;k<SGS.m_bToZ_Axis+1;k++){
				POVDrawUnitCell(pUnit[0], pMap[0], i, j, k, &ar, DSS, ED);
			}

	POVDrawMap(pMap, nbMap, &ar, DSS, LCS);

	ar.WriteString("}\n");
	ar.WriteString("matrix <\n");
	sprintf(szBuffer,"%f , %f , %f ,\n",matrix[0],matrix[1],matrix[2]);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"%f , %f , %f ,\n",matrix[4],matrix[5],matrix[6]);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"%f , %f , %f ,\n",matrix[8],matrix[9],matrix[10]);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"%f , %f , %f  \n",matrix[12],matrix[13],matrix[14]);
	ar.WriteString(szBuffer);
	ar.WriteString(">\n");

	ar.WriteString("scale <1,1,-1>\n");

	ar.WriteString("}\n");

	
	ar.Close();
	f.Close();
	
	
	return true;
}
bool IO::SaveVRMLInput(		CString filename, UnitCell *pUnit[], VoxelMap *pMap[], 
							DisplayStyleStruct DSS, 
							ElementDescription ED[255], 
							SpaceGroupSetup SGS,
							LevelControlSetup LCS, 
							float matrix[16])
{
	char szBuffer[250];
	CString strInputData;
	int nbUnit=0, nbMap=0;

	for(int i=0;i<2;i++){
		if(pUnit[i]!=0) nbUnit++;
		if(pMap[i]!=0) nbMap++;
	}
	if((nbUnit==0)&&(nbMap==0)) return false;
	if((pUnit[0]==0)&&(nbUnit!=0)) return false;
	if((pMap[0]==0)&&(nbMap!=0)) return false;

	CFile f( filename, CFile::modeCreate | CFile::modeWrite);
	CArchive ar(&f,CArchive::store);
	VERIFY(::GetModuleFileName(AfxGetInstanceHandle(), szBuffer, _MAX_PATH));
	CString sPath = (CString)szBuffer;
	sPath = sPath.Left(sPath.ReverseFind('\\'));	
	sPath=sPath+"\\MCE_VRML_head.def";
	CFile f_header(sPath, CFile::modeRead);
	CArchive ar_header(&f_header,CArchive::load);

	ar_header.ReadString(strInputData);
	do
	{
		if(strInputData.Left(13)=="//END OF FILE") break;
		ar.WriteString(strInputData + "\n");
	}while(ar_header.ReadString(strInputData));

	ar_header.Close();
	f_header.Close();

	//define colors

	for(int i=0;i<105;i++){
        sprintf(szBuffer,"DEF  mat_%0.2s Material {\n", ED[i].Name.MakeUpper());
		ar.WriteString(szBuffer);
		sprintf(szBuffer,"diffuseColor %f %f %f\n", ED[i].color[0],ED[i].color[1], ED[i].color[2]);
		ar.WriteString(szBuffer);
		sprintf(szBuffer,"ambientColor 0.2 0.2 0.2\nspecularColor 0.8 0.8 0.8\nshininess 0.2\n}\n");
		ar.WriteString(szBuffer);
	}

	//define bond color;
	sprintf(szBuffer,"DEF  mat_bond Material {\n");
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"diffuseColor %f %f %f\n", DSS.BondsColor[0]/255,DSS.BondsColor[1]/255,DSS.BondsColor[2]/255);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"ambientColor 0.2 0.2 0.2\nspecularColor 0.8 0.8 0.8\nshininess 0.2\n}\n");
	ar.WriteString(szBuffer);

	//define map colors;
	for(int j=0;j<nbMap;j++)
		for(int i=0;i<3;i++){

			sprintf(szBuffer,"DEF  mat_%d%d Material {\n", j,i);
			ar.WriteString(szBuffer);
			sprintf(szBuffer,"diffuseColor %f %f %f\n", pMap[j]->layer_color[i][0],pMap[j]->layer_color[i][1],pMap[j]->layer_color[i][2]);
			ar.WriteString(szBuffer);
			ar.WriteString("ambientColor 0.2 0.2 0.2\nspecularColor 0.8 0.8 0.8\nshininess 0.2\n}\n");

		}

	//background color
	sprintf(szBuffer,"DEF BackgroundColor Info { \nstring \"");
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"%0.2f %0.2f %0.2f\"\n}\n", DSS.BgColor[0], DSS.BgColor[1], DSS.BgColor[2]);
	ar.WriteString(szBuffer);

	//cell color
	sprintf(szBuffer,"DEF  mat_Cell Material {\n");
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"diffuseColor %f %f %f\n", 0.9,0.9,0.9);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"ambientColor 0.2 0.2 0.2\nspecularColor 0.8 0.8 0.8\nshininess 0.2\n}\n");
	ar.WriteString(szBuffer);


	//save transform matrix
	ar.WriteString("MatrixTransform {\nmatrix	");
	sprintf(szBuffer,"%f  %f  %f  %f\n",matrix[0],matrix[1],matrix[2],matrix[3]);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"	%f  %f  %f  %f\n",matrix[4],matrix[5],matrix[6],matrix[7]);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"	%f  %f  %f  %f\n",matrix[8],matrix[9],matrix[10],matrix[11]);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"	%f  %f  %f  %f\n}\n",matrix[12],matrix[13],matrix[14],matrix[15]);
	ar.WriteString(szBuffer);

	for(int i=SGS.m_bFromX_Axis;i<SGS.m_bToX_Axis+1;i++) {
		for(int j=SGS.m_bFromY_Axis;j<SGS.m_bToY_Axis+1;j++) {
			for(int k=SGS.m_bFromZ_Axis;k<SGS.m_bToZ_Axis+1;k++) {
				VRMLDrawUnitCell(pUnit[0], pMap[0], i, j, k, &ar, DSS, ED);
			}
        }
    }

	VRMLDrawMap(pMap, nbMap, &ar, DSS, LCS);
	
	ar.WriteString("}");

	return true;
}
CString IO::MakeJanaAnglRestr(Molecule molecule) 
{
    CString result;
	for (int j=0;j<molecule.arestr.size();j++) {
		// anglefix 109.47 0.01 C1 C2 C3;
        CString tmp;
		tmp.Format("anglefix %.3f 0.01 %s %s %s\n",  molecule.arestr[j].a,  molecule.arestr[j].A1.AtomLabel,  molecule.arestr[j].A2.AtomLabel,  molecule.arestr[j].A3.AtomLabel);  
        result += tmp;
	}
    return result;
}
CString IO::MakeShleXDistRestr(Molecule molecule)
{//DFIX 1.899 0.016 BR1 C1
    
    CString result;
    for (int j=0;j<molecule.brestr.size();j++) {
		// distfix 1.5 0.001 C1 C2;
        CString tmp;
		tmp.Format("DFIX %.3f 0.02 %s %s\n", molecule.brestr[j].d, molecule.brestr[j].A1.AtomLabel,  molecule.brestr[j].A2.AtomLabel);  
        result += tmp;
	}
    return result;
}
CString IO::MakeJanaDistRestr(Molecule molecule) 
{
    CString result;
    for (int j=0;j<molecule.brestr.size();j++) {
		// distfix 1.5 0.001 C1 C2;
        CString tmp;
		tmp.Format("distfix %.3f 0.001 %s %s\n", molecule.brestr[j].d, molecule.brestr[j].A1.AtomLabel,  molecule.brestr[j].A2.AtomLabel);  
        result += tmp;
	}
    return result;
}
CString IO::MakeCrystalsDistRestr(Molecule molecule)
{//DIST  1.421, 0.001 = O(701) TO C(501)
    CString result;
    for (int j=0;j<molecule.brestr.size();j++) {
		// distfix 1.5 0.001 C1 C2;
        CString tmp;
		tmp.Format("DIST %.4f, 0.001 = %s TO %s\n", molecule.brestr[j].d, Atom::getAtomLabelForCrystals(molecule.brestr[j].A1.AtomLabel)
                                                                        , Atom::getAtomLabelForCrystals(molecule.brestr[j].A2.AtomLabel));  
        result += tmp;
	}
    return result;
}
CString IO::MakeCrystalsAnglRestr(Molecule molecule)
{//ANGLES 107.644, 0.01 = O(701) TO C(501) TO C(504)
    CString result;
	for (int j=0;j<molecule.arestr.size();j++) {
		// anglefix 109.47 0.01 C1 C2 C3;
        CString tmp;
		tmp.Format("ANGLES %.4f, 0.01 = %s TO %s TO %s\n",  molecule.arestr[j].a,  
                                                      Atom::getAtomLabelForCrystals(molecule.arestr[j].A1.AtomLabel),  
                                                      Atom::getAtomLabelForCrystals(molecule.arestr[j].A2.AtomLabel),  
                                                      Atom::getAtomLabelForCrystals(molecule.arestr[j].A3.AtomLabel));  
        result += tmp;
	}
    return result;
}
void IO::SaveRestr(CString filename, std::vector<Molecule> molecules)
{
	char szBuffer[250];
	CFile f( filename, CFile::modeCreate | CFile::modeWrite);
	CArchive ar(&f,CArchive::store);    

    for(int i=0;i<molecules.size();i++) {
		ar.WriteString(MakeJanaDistRestr(molecules[i]));
        ar.WriteString(MakeJanaAnglRestr(molecules[i]));
	}
	ar.Close();
	f.Close();
	return;
}
bool IO::SavePeaks(CString filename, std::vector<Atom> peaks, bool names)
{
	char szBuffer[250];
	CFile f( filename, CFile::modeCreate | CFile::modeWrite);
	CArchive ar(&f,CArchive::store);
	
	ar.WriteString("#EDIT\n"); //CRYSTAL command

    for (int i=0;i<peaks.size();i++)
	{
		// CRYSTALS formated atom list
        if(!peaks[i].visible) continue;

        if(names) {
			CString label = peaks[i].AtomLabel;
			if(label.GetLength()==0) {
				label.Format("%d", i+1);
			} else {
				label = "";
				for (int j = 0; j < peaks[i].AtomLabel.GetLength(); j++) {
					if (isdigit(peaks[i].AtomLabel[j])) { 
						label += peaks[i].AtomLabel[j];
					}
				}
			}
            sprintf(szBuffer,"ATOM %s  %s   X=%10.6f  Y=%10.6f  Z=%10.6f \n", peaks[i].AtomType, label,
                                                                            peaks[i].X,
																		    peaks[i].Y,
																		    peaks[i].Z);
        } else {
		    sprintf(szBuffer,"ATOM Q  %d   X=%10.6f  Y=%10.6f  Z=%10.6f \n", i+1, peaks[i].X,
																		    peaks[i].Y,
																		    peaks[i].Z);
        }
		ar.WriteString(szBuffer);
	}

	ar.WriteString("END\n"); //CRYSTAL command

	ar.Close();
	f.Close();
	return true;
}
bool IO::SavePeaks_XYZ(CString filename, std::vector<Atom> peaks)
{
    char szBuffer[250];
	CFile f( filename, CFile::modeCreate | CFile::modeWrite);
	CArchive ar(&f,CArchive::store);	
    for (int i=0;i<peaks.size();i++) {
        if(!peaks[i].visible) continue;
        sprintf(szBuffer,"%s  %s  %10.6f  %10.6f  %10.6f \n", peaks[i].AtomLabel, peaks[i].AtomType, peaks[i].X, peaks[i].Y, peaks[i].Z);        
		ar.WriteString(szBuffer);
	}
    ar.Close();
	f.Close();
    return true;
}
void IO::SaveVoidMap(CString filename, MCEVoidMap void_map, Map_type mapType)
{
    char szBuffer[250];
    std::vector<std::vector<std::vector<float>>> map = void_map.getVoidMap();
    if(map.size()==0) {
        AfxMessageBox("Map was not calculated!");
        return;
    }
	CFile f( filename, CFile::modeCreate | CFile::modeWrite);
	CArchive ar(&f,CArchive::store);

    if(mapType==Map_type::GSAS_GRD_FILE) {
        ar.WriteString("MCE void map\n");
        sprintf(szBuffer,"%10.6f  %10.6f  %10.6f %10.6f  %10.6f  %10.6f\n", void_map.getUnitCell().A-void_map.getXStep(), 
                                                                            void_map.getUnitCell().B-void_map.getYStep(), 
                                                                            void_map.getUnitCell().C-void_map.getZStep(), 
                                                                            void_map.getUnitCell().Alpha,
                                                                            void_map.getUnitCell().Beta,
                                                                            void_map.getUnitCell().Gama); 
        ar.WriteString(szBuffer);
        sprintf(szBuffer,"      %d      %d       %d\n", void_map.getNbXSteps(), void_map.getNbYSteps(), void_map.getNbZSteps()); 
        ar.WriteString(szBuffer);
        for(int i=0;i<map.size();i++) {        
            for(int j=0;j<map[i].size();j++) {
                for(int k=0;k<map[i][j].size();k++) {  
                    sprintf(szBuffer," %10.6f\n", map[i][j][k]);
                    ar.WriteString(szBuffer);
                }
            }
        }
        ar.WriteString("\n");
    } else if(mapType==Map_type::CRYSTALS_FOU_FILE) {

    }
    ar.Close();
	f.Close();
}
bool IO::SaveCif(CString filename, UnitCell *pUnit, std::vector<Atom> &peaks)
{
	if(pUnit==0) return false;

	char szBuffer[250];
	CFile f( filename, CFile::modeCreate | CFile::modeWrite);
	CArchive ar(&f,CArchive::store);
	CString tmp;
	tmp =  "#######################################################################\n";
	tmp += "#\n";
	tmp += "#                 This file was created by MCE\n";
	tmp += "#\n";
	tmp += "#######################################################################\n\n";
	ar.WriteString(tmp);
	ar.WriteString("data_I\n");
	ar.WriteString("_symmetry_space_group_name_H-M   ");
	ar.WriteString("'"+pUnit->getSpaceGroup()+"'\n");
	vector<CString> equivpos = pUnit->getEquivPositions();
	for(int i=0;i<equivpos.size();i++) {
		sprintf(szBuffer,"%d %s\n", i, equivpos[i].GetString());
		ar.WriteString(szBuffer);
	}
	sprintf(szBuffer,"_cell_length_a                   %7.4f\n", pUnit->A);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"_cell_length_b                   %7.4f\n", pUnit->B);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"_cell_length_c                   %7.4f\n", pUnit->C);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"_cell_angle_alpha                %7.4f\n", pUnit->Alpha);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"_cell_angle_beta                %7.4f\n", pUnit->Beta);
	ar.WriteString(szBuffer);
	sprintf(szBuffer,"_cell_angle_gamma               %7.4f\n", pUnit->Gama);
	ar.WriteString(szBuffer);
	//sprintf(szBuffer,"_cell_volume                     %5.3f\n", 1000.0);
	//ar.WriteString(szBuffer);

	ar.WriteString("loop_\n");
	ar.WriteString("_atom_site_label\n");
	ar.WriteString("_atom_site_type_symbol\n");
	ar.WriteString("_atom_site_fract_x\n");
	ar.WriteString("_atom_site_fract_y\n");
	ar.WriteString("_atom_site_fract_z\n");

	vector<Atom> atoms = pUnit->getAssym_part_atoms_fract();
	for(int i=0;i<atoms.size();i++) {
		sprintf(szBuffer,"%s %s %7.4f %7.4f %7.4f\n", atoms[i].AtomLabel.GetString(), atoms[i].AtomType.GetString(), atoms[i].X, atoms[i].Y, atoms[i].Z);
		ar.WriteString(szBuffer);
	}
	for (int i=0;i<peaks.size();i++) {  
        if(!peaks[i].visible) continue;
        sprintf(szBuffer,"%s %s %7.4f %7.4f %7.4f\n", peaks[i].AtomLabel.GetString(), peaks[i].AtomType.GetString(), peaks[i].X, peaks[i].Y, peaks[i].Z);        
		ar.WriteString(szBuffer);
	}	
	
    ar.Close();
	f.Close();
    return true;
}
bool IO::SavePeaks_m40(CString filename, std::vector<Atom> peaks)
{
    char szBuffer[250];
	CFile f( filename, CFile::modeCreate | CFile::modeWrite);
	CArchive ar(&f,CArchive::store);
	for (int i=0;i<peaks.size();i++) {  
        if(!peaks[i].visible) continue;
		int len = peaks[i].AtomLabel.GetLength();
		CString fixed_size_label = peaks[i].AtomLabel;
		for(int j=len;j<10;j++) {
			fixed_size_label += " ";
		}
		CString X, Y, Z;
		X.Format("%.6f", peaks[i].X);
		Y.Format("%.6f", peaks[i].Y);
		Z.Format("%.6f", peaks[i].Z);
		if(X.Find("-")==-1) X = " " + X;
		if(Y.Find("-")==-1) Y = " " + Y;
		if(Z.Find("-")==-1) Z = " " + Z;
		CString final_string;
		final_string = fixed_size_label + "1  1     1.000000" + X + Y + Z + "\n 0.040000 0.000000 0.000000 0.000000 0.000000 0.000000      0000000000\n"; 
        //sprintf(szBuffer,"%s1  1     1.000000 %.6f %.6f %.6f\n 0.040000 0.000000 0.000000 0.000000 0.000000 0.000000      0000000000\n", fixed_size_label.GetString(), peaks[i].X, peaks[i].Y, peaks[i].Z); 
		ar.WriteString(final_string.GetString());
	}	
    ar.Close();
	f.Close();
    return true;
}


//////////////////////////////////////////////////////////
//privates
//////////////////////////////////////////////////////////
float IO::ReadGSASFobFile_interpolate_density(PhaseData *phase, FourierHeader *hdr,
                        FourierData *data,
                        float x, float y, float z) 
{
  // 3D linear interpolation of Fourier density in unit
   // cell space into cartesian space.
   //

  float a, b, c, a0, b0, c0, a1, b1, c1, da, db, dc, I[14];
  int i, j, k;

 ReadGSASFobFile_trans_coords(phase->orth2cryst, &a, &b, &c, x, y, z);
 ReadGSASFobFile_cryst_to_grid(phase, hdr, &i, &j, &k, a, b, c);

  // If outside grid then return zero density
   // else interpolate from surrounding grid points.
   //
  if (i < 0 || j < 0 || k < 0) return 0.0;

  I[0] = data[k].rho[i + j*hdr->nxyzt[0]];
  I[1] = data[k+1].rho[i + j*hdr->nxyzt[0]];
  I[2] = data[k].rho[(i+1) + j*hdr->nxyzt[0]];
  I[3] = data[k+1].rho[(i+1) + j*hdr->nxyzt[0]];
  I[4] = data[k].rho[i + (j+1)*hdr->nxyzt[0]];
  I[5] = data[k+1].rho[i + (j+1)*hdr->nxyzt[0]];
  I[6] = data[k].rho[(i+1) + (j+1)*hdr->nxyzt[0]];
  I[7] = data[k+1].rho[(i+1) + (j+1)*hdr->nxyzt[0]];

 ReadGSASFobFile_grid_to_cryst(phase, hdr, &a0, &b0, &c0, i, j, k);
 ReadGSASFobFile_grid_to_cryst(phase, hdr, &a1, &b1, &c1, i+1, j+1, k+1);

  da = (a - a0)/(a1 - a0);
  db = (b - b0)/(b1 - b0);
  dc = (c - c0)/(c1 - c0);

  I[8] = (I[1] - I[0])*dc + I[0];
  I[9] = (I[3] - I[2])*dc + I[2];
  I[10] = (I[5] - I[4])*dc + I[4];
  I[11] = (I[7] - I[6])*dc + I[6];

  I[12] = (I[10] - I[8])*db + I[8];
  I[13] = (I[11] - I[9])*db + I[9];

  return ((I[13] - I[12])*da + I[12]);
} // End of interpolate_density(). //

void IO::ReadGSASFobFile_trans_coords(float *m,
		              float *x, float *y, float *z,
		              float a, float b, float c) {
  // Multiply a vector (a, b, c) by matrix m to
   // give a transformed vector (x, y, z).
   //

  *x = m[0]*a + m[1]*b + m[2]*c;
  *y = m[3]*a + m[4]*b + m[5]*c;
  *z = m[6]*a + m[7]*b + m[8]*c;
} // End of trans_coords(). //

void IO::ReadGSASFobFile_cryst_to_grid(PhaseData *phase, FourierHeader *hdr,
		               int *i, int *j, int *k,
 		               float a, float b, float c) {

  // Converts fractional unit cell coordinates (a, b, c) 
   // to the grid point nearest the origin that has
   // fractional coordinates closest to (a, b, c). 
   // If the point lies outside grid in any direction then
   // i, j, and k are all set to -1.
   //

  switch (hdr->msect) {
  case 1:
    *i = (int) floor(b*hdr->nxyzi[0] - hdr->nxyzo[0]);
    *j = (int) floor(c*hdr->nxyzi[1] - hdr->nxyzo[1]);
    *k = (int) floor(a*hdr->nxyzi[2] - hdr->nxyzo[2]);
    break;
  case 2:
    *i = (int) floor(c*hdr->nxyzi[0] - hdr->nxyzo[0]);
    *j = (int) floor(a*hdr->nxyzi[1] - hdr->nxyzo[1]);
    *k = (int) floor(b*hdr->nxyzi[2] - hdr->nxyzo[2]);
    break;
  case 3:
    *i = (int) floor(a*hdr->nxyzi[0] - hdr->nxyzo[0]);
    *j = (int) floor(b*hdr->nxyzi[1] - hdr->nxyzo[1]);
    *k = (int) floor(c*hdr->nxyzi[2] - hdr->nxyzo[2]);
    break;
  default:
    printf("Error: unknown section direction.\n");
    exit(1);
    break;
  }

  if (*i < 0 || *j < 0 || *k < 0 ||
      *i >= (hdr->nxyzt[0] - 1) ||
      *j >= (hdr->nxyzt[1] - 1) || 
      *k >= (hdr->nxyzt[2] - 1)) {
    *i = -1;
    *j = -1;
    *k = -1;
  }
} // End of cryst_to_grid(). //

void IO::ReadGSASFobFile_grid_to_cryst(PhaseData *phase, FourierHeader *hdr,
		               float *a, float *b, float *c,
		               int i, int j, int k) {

  // Convert grid point coordinates (i, j, k) into  
   // fractional coordinates (a, b, c) w.r.t the unit cell
   //axes.
   //

  switch (hdr->msect) {
  case 1:
    *a = (float) (k + hdr->nxyzo[2])/hdr->nxyzi[2];
    *b = (float) (i + hdr->nxyzo[0])/hdr->nxyzi[0];
    *c = (float) (j + hdr->nxyzo[1])/hdr->nxyzi[1];
    break;
  case 2:
    *a = (float) (j + hdr->nxyzo[1])/hdr->nxyzi[1];
    *b = (float) (k + hdr->nxyzo[2])/hdr->nxyzi[2];
    *c = (float) (i + hdr->nxyzo[0])/hdr->nxyzi[0];
    break;
  case 3:
    *a = (float) (i + hdr->nxyzo[0])/hdr->nxyzi[0];
    *b = (float) (j + hdr->nxyzo[1])/hdr->nxyzi[1];
    *c = (float) (k + hdr->nxyzo[2])/hdr->nxyzi[2];
    break;
  default:
    printf("Error: unknown section direction.\n");
    exit(1);
    break;
  }
} // End of grid_to_cryst(). //


void IO::ReadGSASFobFile_output_crystals(PhaseData *phase, FourierHeader *hdr,FourierData *foudata, int interactive, UnitCell *pUnit, VoxelMap *pMap)
{
	#define PI 3.1415926536
	int i, j, k, nx, ny, nz, nxy;
	float x, y, z;
	float xstart, ystart, zstart;
	float xstep, ystep, zstep;
	float xend, yend, zend;
	
	CString strInputData;

	pMap->mType=GSAS_FOB_FILE; // Map type identification

	ReadGSASFobFile_setup_trans_matrices(phase);

	ReadGSASFobFile_get_box_limits(	phase, hdr,
									&xstart, &xstep, &xend,
									&ystart, &ystep, &yend,
									&zstart, &zstep, &zend);
	int loop=0;
    float m[3][3];
	for (i=0;i<3;i++)
	{
		for(j=0;j<3;j++)
		{
			m[i][j] = phase->cryst2orth[loop];
			loop++;
		}
	}
    pUnit->setFtC_Matrix(m);
	for (i=0;i<3;i++)
	{
		pUnit->Translation_Matrix[i]=phase->cryst2orth[9+i];
	}

	pUnit->A = phase->cell[0];
	pUnit->B = phase->cell[1];
	pUnit->C = phase->cell[2];
	pUnit->Alpha = phase->cell[3];
	pUnit->Beta = phase->cell[4];
	pUnit->Gama = phase->cell[5];

	pMap->AxisParam[0].Start=xstart;
	pMap->AxisParam[0].Step=xstep;
	pMap->AxisParam[0].End=xend;
	pMap->AxisParam[0].Division=1.0;
    pMap->AxisParam[0].Dimension = phase->cell[0];

	pMap->AxisParam[1].Start=ystart;
	pMap->AxisParam[1].Step=ystep;
	pMap->AxisParam[1].End=yend;
	pMap->AxisParam[1].Division=1.0;
    pMap->AxisParam[1].Dimension = phase->cell[1];

	pMap->AxisParam[2].Start=zstart;
	pMap->AxisParam[2].Step=zstep;
	pMap->AxisParam[2].End=zend;
	pMap->AxisParam[2].Division=1.0;
    pMap->AxisParam[2].Dimension = phase->cell[2];

	nx = (int) ceil((xend - xstart)/xstep);
	ny = (int) ceil((yend - ystart)/ystep);
	nz = (int) ceil((zend - zstart)/zstep);
  
	pMap->Lines=nx;
	pMap->Columns=ny;
	pMap->Sections=nz;

	if ((pMap->Lines>MAX_OPTIMAL_VOXEL_DIEMENSION)|
		(pMap->Columns>MAX_OPTIMAL_VOXEL_DIEMENSION)|
		(pMap->Sections>MAX_OPTIMAL_VOXEL_DIEMENSION)){
		AfxMessageBox("Too large ELD field(in one direction bigger than 400)!\n Program will use more memory...\n");
	}

	pMap->MemoryForVoxel();

	nxy = nx*ny;
	for (k = 0, z = zstart; k < nz; k++, z += zstep) {
		for (i = 0, x = xstart; i < nx; i++, x += xstep) {
			  for (j = 0, y = ystart; j < ny; j++, y += ystep) {
				pMap->Voxel[i][j][k]=ReadGSASFobFile_interpolate_density(phase, hdr, foudata, x, y, z);
			  }
		}
	}

	if(phase->natoms>IO_MAXATOMS) {
		AfxMessageBox("Error: There are more than 2000 atoms in input file.");
		return;
	}

	if (phase->natoms>0)
	{
		int atom_record_lenght=14;
		for (i=0;i<phase->natoms;i++)
		{
			CString tmpName = phase->type[i];
            pUnit->AddNewAtomInFract(Atom::getAtomType(tmpName), tmpName, phase->x[i], phase->y[i], phase->z[i]);
		}
	}

}
void IO::ReadGSASFobFile_get_box_limits(PhaseData *phase, FourierHeader *hdr,
		                float *xstart, float *xstep, float *xend,
		                float *ystart, float *ystep, float *yend,
		                float *zstart, float *zstep, float *zend) {

  // Get the minimum and maximum extent in cartesian coordinates
   //of the GSAS grid, which is in unit cell coordinates.
   //
   
  int ci, cj, ck, i, j, k;
  float a, b, c, x, y, z;
  float xmin, ymin, zmin, xmax, ymax, zmax;

  xmin = 10000000.0;
  ymin = 10000000.0;
  zmin = 10000000.0;
  xmax = -10000000.0;
  ymax = -10000000.0;
  zmax = -10000000.0;
  
  for (ck = 0; ck <= 1; ck++) {
    for (cj = 0; cj <= 1; cj++) {
      for (ci = 0; ci <= 1; ci++) {
	      i = ci*(hdr->nxyzt[0] - 1);
 	      j = cj*(hdr->nxyzt[1] - 1);
	      k = ck*(hdr->nxyzt[2] - 1);
	     ReadGSASFobFile_grid_to_cryst(phase, hdr, &a, &b, &c, i, j, k);
	     ReadGSASFobFile_trans_coords(phase->cryst2orth, &x, &y, &z, a, b, c);
	      xmin = (x < xmin) ? x : xmin;
	      xmax = (x > xmax) ? x : xmax;
	      ymin = (y < ymin) ? y : ymin;
	      ymax = (y > ymax) ? y : ymax;
	      zmin = (z < zmin) ? z : zmin;
	      zmax = (z > zmax) ? z : zmax;
      }
    }
  }

  *xstart = xmin;
  *xstep = phase->cell[0]/hdr->nxyzi[0];
  *xend = xmax;
  *ystart = ymin;
  *ystep = phase->cell[1]/hdr->nxyzi[1];
  *yend = ymax;
  *zstart = zmin;
  *zstep = phase->cell[2]/hdr->nxyzi[2];
  *zend = zmax;
} // End of get_box_limits(). //


void IO::ReadGSASFobFile_setup_trans_matrices(PhaseData *phase) 
{

  // m[0] == m11, m[1] = m12, m[2] = m13
   // m[3] == m21, m[4] = m22, m[5] = m23
   // m[6] == m31, m[7] = m32, m[8] = m33
   // m[9] == x trans, m[10] = y trans, m[11] = z trans
   // phase->cryst2orth and phase->orth2cryst have the
   // layout above.
   // phase->cryst2orth is the matrix with which to multiply 
   // fractional unit cell coordinates to get coordinates in
   // an orthonormal coordinate system with x along a,
   // z along c* and y normal to the x-z plane.
   // phase->orth2cryst is the inverse matrix to do the
   // inverse operation.
   //
	float a, b, c, alpha, beta, gamma;
	float cos_alpha, sin_alpha, cos_beta, sin_beta, cos_gamma, sin_gamma, vol;
    float a_star, b_star, c_star, cos_alpha_star, cos_beta_star;

  a = phase->cell[0];
  b = phase->cell[1];
  c = phase->cell[2];
  alpha = phase->cell[3];
  beta = phase->cell[4];
  gamma = phase->cell[5];

  cos_alpha = (float) cos(alpha*PI/180.0);
  sin_alpha = (float) sin(alpha*PI/180.0);
  cos_beta = (float) cos(beta*PI/180.0);
  sin_beta = (float) sin(beta*PI/180.0);
  cos_gamma = (float) cos(gamma*PI/180.0);
  sin_gamma = (float) sin(gamma*PI/180.0);

  vol = a*a*b*b*c*c*
        (1 - cos_alpha*cos_alpha - cos_beta*cos_beta - cos_gamma*cos_gamma
	 + 2*cos_alpha*cos_beta*cos_gamma);
  if (vol <= 0.0) {
    printf("Error: Bad unit cell volume: vol^2 = %f\n", vol);
    printf("Unit cell params have not been read correctly.\n");
    printf("You can still output in ASCII format.\n");
    exit(1);
  }

  vol = sqrt(vol);

  a_star = b*c*sin_alpha/vol;
  b_star = a*c*sin_beta/vol;
  c_star = a*b*sin_gamma/vol;
  cos_alpha_star = (cos_beta*cos_gamma - cos_alpha)/(sin_beta*sin_gamma);
  cos_beta_star = (cos_alpha*cos_gamma - cos_beta)/(sin_alpha*sin_gamma);

  phase->cryst2orth[0] = (float) a;
  phase->cryst2orth[1] = (float) b*cos_gamma;
  phase->cryst2orth[2] = (float) c*cos_beta;
  phase->cryst2orth[3] = (float) 0.0;
  phase->cryst2orth[4] = (float) b*sin_gamma;
  phase->cryst2orth[5] = (float) -c*sin_beta*cos_alpha_star;
  phase->cryst2orth[6] = (float) 0.0;
  phase->cryst2orth[7] = (float) 0.0;
  phase->cryst2orth[8] = (float) 1.0/c_star;
  phase->cryst2orth[9] = 0.0;
  phase->cryst2orth[10] = 0.0;
  phase->cryst2orth[11] = 0.0;

  phase->orth2cryst[0] = (float) 1.0/a;
  phase->orth2cryst[1] = (float) -cos_gamma/(a*sin_gamma);
  phase->orth2cryst[2] = (float) a_star*cos_beta_star;
  phase->orth2cryst[3] = (float) 0.0;
  phase->orth2cryst[4] = (float) 1.0/(b*sin_gamma);
  phase->orth2cryst[5] = (float) b_star*cos_alpha_star;
  phase->orth2cryst[6] = (float) 0.0;
  phase->orth2cryst[7] = (float) 0.0;
  phase->orth2cryst[8] = (float) c_star;
  phase->orth2cryst[9] = (float) 0.0;
  phase->orth2cryst[10] = (float) 0.0;
  phase->orth2cryst[11] = (float) 0.0;
} // End of setup_trans_matrices(). //


void IO::ReadGSASFobFile_get_gsas_line(char *line, FILE *file) {

  char ch;
  int nRead;
    
  //The line parameter must point to a character array
   //of at least 81 bytes, enough for 80 characters 
   // followed by a NULL.
   //
   //
   // Routine that reads lines from all types of GSAS 
   // ASCII files (UNIX and PC). 
   //
  nRead = fread(line, 1, 80, file);
  if (nRead == 80) {
    fread(&ch, 1, 1, file);
    switch (ch) {
    case '\n':
      // Don't do anything. //
      break;
    case '\r':
      // Jump forward one char because this is 
       // probably a PC format text file with \r\n
       // pairs terminating the lines. 
       //
      fseek(file, 1, SEEK_CUR);
      break;
    default:
      //If there is no \n or \r after 80 chars then
       //this is probably a UNIX format GSAS file.
       //
      fseek(file, -1, SEEK_CUR);
      break;
    }
  }

  line[nRead] = '\0';
} // End of get_gsas_line(). //

void IO::ReadGSASFobFile_parse_expfile(CString filename, PhaseData *data, int interactive) 
{
  FILE *exp_file;
  char line[81], exp_filename[255], dummy_string[80];
  int i, nphases, phase, natoms[10];
  char phase_name[10][80];
  float cell[10][6];
  float xyzo[10][1000][4];
  char type[10][1000][3];
  char label[10][1000][12];
  int  ch;

	ch=filename.GetLength();
    filename.SetAt( ch-3, 'e');
	filename.SetAt( ch-2, 'x');
	filename.SetAt( ch-1, 'p');
  
  // Try to open the *.EXP file. //

  if ((exp_file = fopen(filename, "rb")) == NULL) {
    
    // Can't open with lower case filename so try all upper case. // 
    filename.SetAt( ch-3, 'E');
	filename.SetAt( ch-2, 'X');
	filename.SetAt( ch-1, 'P');
  
    if ((exp_file = fopen(exp_filename, "rb")) == NULL) {
	  MessageBox(NULL, "Can't find *.exp or *.EXP file", "WinMain", MB_ICONEXCLAMATION); 
      data->cell[0] = data->cell[1] = data->cell[2] = 0.0;
      data->cell[2] = data->cell[3] = data->cell[4] = 0.0;
      data->natoms = 0;
      return;
    }
  }
  
  // Parse the *.EXP file to get the phase information. //
  nphases = 0;
ReadGSASFobFile_get_gsas_line(line, exp_file);
  while (!feof(exp_file) && nphases < 10) {
    if ((strncmp(&line[4], "    PNAM", 8)) == 0) {
      // Read name. //
      sscanf(&line[14], "%s", phase_name[nphases]);

      // Read number of atoms. //
    ReadGSASFobFile_get_gsas_line(line, exp_file);
      while (strncmp(&line[6], " NATOM", 6) != 0) {
	     ReadGSASFobFile_get_gsas_line(line, exp_file);
      }
      sscanf(&line[14], "%d", &natoms[nphases]);
      if (natoms[nphases] > 1000) {
        fprintf(stderr,
		"Error: Phase %d has %d atoms.\n", nphases, natoms[nphases]);
        fprintf(stderr, "Can't handle more than %d atoms.\n", IO_MAXATOMS);
        fprintf(stderr, "Recompile after increasing MAXATOMS.\n");
        exit(1);
      }

      // Read a, b, c. //
     ReadGSASFobFile_get_gsas_line(line, exp_file);
      while (strncmp(&line[6], "ABC", 2) != 0) {
	    ReadGSASFobFile_get_gsas_line(line, exp_file);
      }
      sscanf(&line[12], "%f%f%f%s",
	           &cell[nphases][0], &cell[nphases][1], &cell[nphases][2],
	           dummy_string);

      // Read alpha, beta, gamma. //
     ReadGSASFobFile_get_gsas_line(line, exp_file);
      while (strncmp(&line[6], "ANGLES", 6) != 0) {
	     ReadGSASFobFile_get_gsas_line(line, exp_file);
      }
      sscanf(&line[12], "%f%f%f",
	           &cell[nphases][3], &cell[nphases][4], &cell[nphases][5]);

      // Read atom info. //
     ReadGSASFobFile_get_gsas_line(line, exp_file);
      while (strncmp(&line[6], "AT", 2) != 0) {
	     ReadGSASFobFile_get_gsas_line(line, exp_file);
      }
      for (i = 0; i < natoms[nphases]; i++) {
	strncpy(type[nphases][i], &line[14], 2);
	type[nphases][i][2] = '\0';
        strncpy(label[nphases][i], &line[62], 11);
        label[nphases][i][11] = '\0';
        line[62] = '\0';
        sscanf(&line[22], "%f%f%f%f", 
               &xyzo[nphases][i][0], &xyzo[nphases][i][1],
               &xyzo[nphases][i][2], &xyzo[nphases][i][3]);
       ReadGSASFobFile_get_gsas_line(line, exp_file);
	     ReadGSASFobFile_get_gsas_line(line, exp_file);
      }
      nphases++;
    }
   ReadGSASFobFile_get_gsas_line(line, exp_file);
  }
  
  if (interactive == 0) 
    printf("Number of phases found in EXP file: %d\n", nphases);

  if (nphases == 0) {
    if (interactive == 0) {
      printf("No phase information. Setting unit cell to zero.\n");
      printf("This may make sone output formats unreadable.\n");
    }
    data->cell[0] = data->cell[1] = data->cell[2] = 0.0;
    data->cell[2] = data->cell[3] = data->cell[4] = 0.0;
    data->natoms = 0;
    return;
  }

  phase = 0;
  if (interactive == 0) {
    for (i = 0; i < nphases; i++) {
      printf("Phase %d: Phase Name = %s\n", (i+1), phase_name[i]);
      printf("Cell: %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f\n",
	     cell[i][0], cell[i][1], cell[i][2], 
	     cell[i][3], cell[i][4], cell[i][5]);
      printf("Number of atoms in asymmetric unit: %d\n", natoms[i]);
    }
  
    if (nphases > 1) {
      do {
	printf("Which phase is the fourier map from (1 - %d): ", nphases);
	scanf("%d", &phase);
	phase--;
      } while(phase <= 0 || phase > nphases);
    } else {
      phase = 0;
    }
  }

  data->cell[0] = cell[phase][0];
  data->cell[1] = cell[phase][1];
  data->cell[2] = cell[phase][2];
  data->cell[3] = cell[phase][3];
  data->cell[4] = cell[phase][4];
  data->cell[5] = cell[phase][5];
  data->natoms = natoms[phase];
  for (i = 0; i < natoms[phase]; i++) {
    data->x[i] = xyzo[phase][i][0];
    data->y[i] = xyzo[phase][i][1];
    data->z[i] = xyzo[phase][i][2];
    data->occ[i] = xyzo[phase][i][3];
    strcpy(data->type[i], type[phase][i]);
    strcpy(data->label[i], label[phase][i]);
  }

  fclose(exp_file);
  return;
 // End of parse_expfile(). //

}
bool IO::ReadXDGRDFile(CString filename, UnitCell *pUnit, VoxelMap *pMap)
{	
	float lines, columns, sections;
	float start_x, start_y, start_z;
	float end_x, end_y, end_z;
	int i, j, k, atoms;
	char szString[10], t[10];
	float x, y, z;

	CString strInputData;
	
	CFile f(filename, CFile::modeRead);
	CArchive ar(&f,CArchive::load);
	
	do{
		if(!ar.ReadString(strInputData)) return false;//skip header
	}while((strInputData!="! Gridpoints, Origin, Physical Dimensions"));
	
	ar.ReadString(strInputData); // number of section  
	sscanf(strInputData,"%f	%f	%f", &lines, &columns, &sections);

	pMap->Lines = (int) lines;
	pMap->Columns = (int) columns;
	pMap->Sections = (int) sections;

	if ((pMap->Lines>MAX_OPTIMAL_VOXEL_DIEMENSION)|
		(pMap->Columns>MAX_OPTIMAL_VOXEL_DIEMENSION)|
		(pMap->Sections>MAX_OPTIMAL_VOXEL_DIEMENSION)){
		AfxMessageBox("Too large ELD field(in one direction bigger than 400)!\n Program will use more memory...\n");
	}
	
	pMap->MemoryForVoxel();

	ar.ReadString(strInputData);//start
	sscanf(strInputData,"%f %f %f",&start_x, &start_y, &start_z);
	
	ar.ReadString(strInputData);//end
	sscanf(strInputData,"%f %f %f", &end_x, &end_y, &end_z);
	
	pMap->AxisParam[0].Start = start_x;
	pMap->AxisParam[0].Step = fabs(end_x - start_x)/(pMap->Lines - 1);
	pMap->AxisParam[0].End = end_x;
	pMap->AxisParam[0].Division = 1;
    pMap->AxisParam[0].Dimension = end_x;
	
	pMap->AxisParam[1].Start = start_y;
	pMap->AxisParam[1].Step = fabs(end_y - start_y)/(pMap->Columns - 1);
	pMap->AxisParam[1].End = end_y;
	pMap->AxisParam[1].Division = 1;
    pMap->AxisParam[1].Dimension = end_y;
	
	if (pMap->mType == XD_GRD2D_FILE){
		pMap->AxisParam[2].Start = start_z;
		pMap->AxisParam[2].Step = 0;
		pMap->AxisParam[2].End = end_z;
		pMap->AxisParam[2].Division = 1;
        pMap->AxisParam[2].Dimension = end_z;
	}else{
		pMap->AxisParam[2].Start = start_z;
		pMap->AxisParam[2].Step = fabs(end_z - start_z)/(pMap->Sections - 1);
		pMap->AxisParam[2].End = end_z;
		pMap->AxisParam[2].Division = 1;
        pMap->AxisParam[2].Dimension = end_z;
	}
	
	do{
		if(!ar.ReadString(strInputData)) return false;
	}while((strInputData!="! Objects"));

	ar.ReadString(strInputData);// 
	sscanf(strInputData,"%d", &atoms);
	

	if(atoms>IO_MAXATOMS) {
		AfxMessageBox("Error: There are more than 2000 atoms in input file.");
		return false;
	}

	for (i=0;i<atoms;i++){
		ar.ReadString(strInputData);// reading atoms
		sscanf(strInputData,"%s %f %f %f %s", &szString, &x, &y, &z, &t);
		
			//strncpy(pMolecule->a[i].name,szString,2);        
		pUnit->AddNewAtomInFract(Atom::getAtomType(szString), szString, x, y, z);		
	}
	TRACE("Finished reading Molecule data");

	do{
		if(!ar.ReadString(strInputData)) return false;
	}while((strInputData!="! Values"));
	
	if (pMap->mType == XD_GRD2D_FILE)
	{
		for (j=0;j<pMap->Columns;j++)
		{
			for (i=0;i<pMap->Lines;i++)
			{
				pMap->Voxel[i][j][0] = GetNumberFromFile(&ar);
			}
		}
	}
	if (pMap->mType == XD_GRD3D_FILE)
	{
		for (k=0;k<pMap->Sections;k++)
		{
			for (j=0;j<pMap->Columns;j++)
			{
				for (i=0;i<pMap->Lines;i++)
				{
					pMap->Voxel[i][j][k] = GetNumberFromFile(&ar);
				}
			}
		}
	}
	TRACE("Finished reading Voxel Map");
	ar.Close();
	f.Close();

	pUnit->SetSpaceGroupSpecial("P1");
	pUnit->Option_dummy_fill = true;
	pUnit->Option_setToUnitCell = false;
	return true;

}
float IO::GetNumberFromFile(CArchive *ar)
{
	char tmp[1];
	CString number="";
	float result;
    //seek begin
	do{
		if(ar->Read(tmp, 1)<1) {
            return 0;            
        }        
	}while(((tmp[0]==32)|(tmp[0]==13)|(tmp[0]==10)));

	do{
		number += tmp[0];
		if(ar->Read(tmp, 1)<1) break;
	}while(!((tmp[0]==32)|(tmp[0]==13)|(tmp[0]==10)));

	
	sscanf(number, "%f", &result);

	return result;
}
bool IO::GetNumberFromString(CString &str, int nb_char, float &number)
{
    if(str.IsEmpty()) return false;

    CString res = str.Mid(0, 12);
    sscanf(res, "%f", &number);
    str = str.Mid(12, str.GetLength());
    return true;
}
bool IO::ReadLineFromString(CString &str, CString &line)
{
    str.Trim();
    if(str.GetLength() == 0) return false;

    int i = str.Find("\n");
    if(i!=-1) {
        line = str.Mid(0, i);
        str = str.Mid(i+1);
    } else  {
        line = str;
        str = "";
        return true;
    }
    return true;
}
bool IO::GetDefinedVariables(CString str, map<CString, CString> &variables)
{
    //#define MCE_VERSION "3.6.04"
    //#define MCE_DATE "2017/07/21"
    CString line;
    while (ReadLineFromString(str, line)) {
        CString param, value;
        line.Replace("#define ", "");
        GetFirstParamFromString(line, param, value);
        GetStringBetweenQuotes(value, value);
        variables[param] = value;
    }

    return true;
}
bool IO::IsValidInt(const CString& text, long& value)
{
    LPCTSTR ptr = (LPCTSTR) text;
    LPTSTR endptr;
    value = _tcstol(ptr, &endptr, 10);
    return (*ptr && endptr - ptr == text.GetLength());
}

bool IO::IsValidFloat(const CString& text, float& value)
{
    LPCTSTR ptr = (LPCTSTR) text;
    LPTSTR endptr;
    value = _tcstod(ptr, &endptr);
    return (*ptr && endptr - ptr == text.GetLength());
}
bool IO::GetStringBetweenQuotes(CString str, CString &res)
{
    int i = str.Find("\"", 0);
    if(i!=-1) {
        int j = str.Find("\"", i+1);
        if(j!=-1) {
            res = str.Mid(i+1, j-i-1);
        } else {
            return false;
        }
    } else {
        return false;
    }
    return true;
}
bool IO::ReadNormalGRDFile(CString filename, UnitCell *pUnit, VoxelMap *pMap)
{
	CString strInputData;

	// support variable
	int nLineNumber=0;
	int i, j;
	int x,y,z;
	
	TRACE("Reading grd file: %s\n",filename);
	
	CFile f( filename, CFile::modeRead);
	CArchive ar(&f,CArchive::load);

	ar.ReadString(strInputData); // GRD information string
	ar.ReadString(strInputData); // unit cell parameters

	float a,b,c,al,be,ga;
	sscanf(strInputData,"%f %f %f %f %f %f",&a,&b,&c,&al,&be,&ga);
	pUnit->A = a;
	pUnit->B = b;
	pUnit->C = c;
	pUnit->Alpha = al;
	pUnit->Beta = be;
	pUnit->Gama = ga;    

	ar.ReadString(strInputData); // steps in lattice directions
	int stepsx,stepsy,stepsz;
	sscanf(strInputData,"%d %d %d",&stepsx,&stepsy,&stepsz);

    /*
    if(pMap->isPeriodic) {
        stepsx++;
        stepsy++;
        stepsz++;
    }
    */

	// let simulate start at 0 and 0.1 strps

	for (i=0;i<3;i++)
	{
		pMap->AxisParam[i].Start=0;		
		pMap->AxisParam[i].End=1;
		pMap->AxisParam[i].Division=stepsx;
		
	}
	pMap->AxisParam[0].Step=a/((float)stepsx-1);
	pMap->AxisParam[1].Step=b/((float)stepsy-1);
	pMap->AxisParam[2].Step=c/((float)stepsz-1);
   	pMap->AxisParam[0].Dimension=a;
	pMap->AxisParam[1].Dimension=b;
	pMap->AxisParam[2].Dimension=c;


	pMap->Lines=stepsx;
	pMap->Columns=stepsy;
	pMap->Sections=stepsz;

	if ((pMap->Lines>MAX_OPTIMAL_VOXEL_DIEMENSION)|
		(pMap->Columns>MAX_OPTIMAL_VOXEL_DIEMENSION)|
		(pMap->Sections>MAX_OPTIMAL_VOXEL_DIEMENSION)){
		AfxMessageBox("Too large ELD field(in one direction bigger than 400)!\n Program will use more memory...\n");
	}

	pMap->MemoryForVoxel();
/*
    if(pMap->isPeriodic) {
        for(x=0;x<pMap->Lines;x++)
	    {
		    for(y=0;y<pMap->Columns;y++)
		    {
			    for (z=0;z<pMap->Sections;z++)
			    {
                    if((x==(pMap->Lines-1)) || (y==(pMap->Columns-1)) || (z==(pMap->Sections-1))) {
                        int i=x, j=y, k=z;
                        if(x == (pMap->Lines-1)) {
		                    i=0;
	                    }
	                    if(y == (pMap->Columns-1)) {
		                    j=0;
	                    }
                        if(z == (pMap->Sections-1)) {
		                    k=0;
	                    }
                        pMap->Voxel[x][y][z] = pMap->Voxel[i][j][k];
                    } else {
				        pMap->Voxel[x][y][z] = GetNumberFromFile(&ar);
                    }
			    }
		    }
	    }
        
    } else {
    */
	    for(x=0;x<pMap->Lines;x++) {
		    for(y=0;y<pMap->Columns;y++) {
			    for (z=0;z<pMap->Sections;z++) {
				    pMap->Voxel[x][y][z] = GetNumberFromFile(&ar);
			    }
		    }
	    }
    //}

	TRACE("Finished reading Voxel Map");
  
	ar.Close();
	f.Close();

	pUnit->SetSpaceGroup("P1");
    float m1[3][3], m2[3][3];
    pUnit->getFTC_Matrix(m1);
    pUnit->getCTF_Matrix(m2);
	for(i=0;i<3;i++)
		for(j=0;j<3;j++){
            pMap->MatrixFtC[i][j] = m1[i][j];
            pMap->MatrixCtF[i][j] = m2[i][j];
		}
	return true;
}
void IO::POVDrawUnitCell(UnitCell *pUnit, VoxelMap *pMap, int x, int y, int z, CArchive *ar, DisplayStyleStruct DSS, ElementDescription ED[255])
{
	float Trans_xyz[3];
	char szBuffer[400];

	//compute real vector from x, y, z
	Trans_xyz[0] = (float) x * (pUnit->GetUnitPoint_cart(1).X - pUnit->GetUnitPoint_cart(0).X)
				 + (float) y * (pUnit->GetUnitPoint_cart(2).X - pUnit->GetUnitPoint_cart(0).X)
				 + (float) z * (pUnit->GetUnitPoint_cart(4).X - pUnit->GetUnitPoint_cart(0).X);
	Trans_xyz[1] = (float) x * (pUnit->GetUnitPoint_cart(1).Y - pUnit->GetUnitPoint_cart(0).Y) 
				 + (float) y * (pUnit->GetUnitPoint_cart(2).Y - pUnit->GetUnitPoint_cart(0).Y)
				 + (float) z * (pUnit->GetUnitPoint_cart(4).Y - pUnit->GetUnitPoint_cart(0).Y);
	Trans_xyz[2] = (float) x * (pUnit->GetUnitPoint_cart(1).Z - pUnit->GetUnitPoint_cart(0).Z) 
				 + (float) y * (pUnit->GetUnitPoint_cart(2).Z - pUnit->GetUnitPoint_cart(0).Z) 
				 + (float) z * (pUnit->GetUnitPoint_cart(4).Z - pUnit->GetUnitPoint_cart(0).Z);

	if(pMap!=0)
	if( (pMap->mType == CRYSTALS_FOU_FILE)|(pMap->mType == GSAS_GRD_FILE)){

		//pUnit->SetAllBasicMoleculesInCart();
		if(!pUnit->FillInUnitCell()){
			AfxMessageBox("FillInUnitCell returns false");
			return;
		}
	}else{
		if(!pUnit->FillInUnitCell()){
			AfxMessageBox("FillInUnitCell returns false");
			return;
		}
	}
	
	//output atoms
	if((DSS.DisplayStyle!=Off)&&(DSS.DisplayStyle!=Line))
		if(!pUnit->Option_Organic){
            std::vector<Atom> UCAtoms = pUnit->getUnitcell_atoms_cart();
            for(int i=0;i<UCAtoms.size();i++){

				sprintf(szBuffer,"sphere {< %f , %f , %f > , ",	UCAtoms[i].X + Trans_xyz[0],
																UCAtoms[i].Y + Trans_xyz[1],
																UCAtoms[i].Z + Trans_xyz[2]);
				ar->WriteString(szBuffer);	

				if ((DSS.DisplayStyle==Ball_and_line)|
					(DSS.DisplayStyle==Ball_and_cylinder)){
					ar->WriteString("AtomRadius");
				}
				if(DSS.DisplayStyle==Cylinder){
					ar->WriteString("bondcylRadius");
				}
				if (DSS.DisplayStyle==Scaled_ball_and_cylinder){
					sprintf(szBuffer,"AtomRadiusFor%0.2s*ScaledFactor",UCAtoms[i].AtomType.MakeUpper());
					ar->WriteString(szBuffer);
				}

				sprintf(szBuffer," pigment {color ColorFor%0.2s}}\n",UCAtoms[i].AtomType.MakeUpper());
				ar->WriteString(szBuffer);
			}
		}else{
            std::vector<Molecule> UCMolecule = pUnit->getUnitcell_molecule_cart();
            for(int j=0;j<UCMolecule.size();j++){
				for(int i=0;i<UCMolecule[j].GetAtomNb();i++){

					sprintf(szBuffer,"sphere {< %f , %f , %f > , ",	UCMolecule[j].GetAtom(i).X + Trans_xyz[0],
																	UCMolecule[j].GetAtom(i).Y + Trans_xyz[1],
																	UCMolecule[j].GetAtom(i).Z + Trans_xyz[2]);
					ar->WriteString(szBuffer);	

					if ((DSS.DisplayStyle==Ball_and_line)|
						(DSS.DisplayStyle==Ball_and_cylinder)){
						ar->WriteString("AtomRadius");
					}
					if(DSS.DisplayStyle==Cylinder){
						ar->WriteString("bondcylRadius");
					}
					if (DSS.DisplayStyle==Scaled_ball_and_cylinder){
						sprintf(szBuffer,"AtomRadiusFor%0.2s*ScaledFactor",UCMolecule[j].GetAtom(i).AtomType.MakeUpper());
						ar->WriteString(szBuffer);
					}

					sprintf(szBuffer," pigment {color ColorFor%0.2s}}\n",UCMolecule[j].GetAtom(i).AtomType.MakeUpper());
					ar->WriteString(szBuffer);
				}
			}
		}
	//output bonds
	if((DSS.DisplayStyle!=Off)&&(DSS.bBondsVisible))
		if(!pUnit->Option_Organic){
			for(int i=0;i<pUnit->Bond.size();i++){
                if(UnitCell::Distance(pUnit->Bond[i].A1X + Trans_xyz[0],
                                      pUnit->Bond[i].A1Y + Trans_xyz[1],
                                      pUnit->Bond[i].A1Z + Trans_xyz[2],
                                      pUnit->Bond[i].A2X + Trans_xyz[0],
                                      pUnit->Bond[i].A2Y + Trans_xyz[1],
                                      pUnit->Bond[i].A2Z + Trans_xyz[2]) > 0.01) {
				    if(DSS.iBondsColorStyle==TwoColors)
				    {
					    sprintf(szBuffer,"cylinder{< %f, %f, %f > , < %f, %f, %f>, ",	pUnit->Bond[i].A1X + Trans_xyz[0],
																					    pUnit->Bond[i].A1Y + Trans_xyz[1],
																					    pUnit->Bond[i].A1Z + Trans_xyz[2], 
																					    ((pUnit->Bond[i].A1X+pUnit->Bond[i].A2X)/2) + Trans_xyz[0],
																					    ((pUnit->Bond[i].A1Y+pUnit->Bond[i].A2Y)/2) + Trans_xyz[1],
																					    ((pUnit->Bond[i].A1Z+pUnit->Bond[i].A2Z)/2) + Trans_xyz[2]);
					    ar->WriteString(szBuffer);
					    ar->WriteString("bondcylRadius");
					    sprintf(szBuffer," pigment {color ColorFor%0.2s}}\n",ED[(int)pUnit->Bond[i].A1AtNb-1].Name);
					    ar->WriteString(szBuffer);

					    sprintf(szBuffer,"cylinder{< %f, %f, %f > , < %f, %f, %f>, bondcylRadius",	pUnit->Bond[i].A2X + Trans_xyz[0],
																								    pUnit->Bond[i].A2Y + Trans_xyz[1],
																								    pUnit->Bond[i].A2Z + Trans_xyz[2], 
																								    ((pUnit->Bond[i].A1X+pUnit->Bond[i].A2X)/2) + Trans_xyz[0],
																								    ((pUnit->Bond[i].A1Y+pUnit->Bond[i].A2Y)/2) + Trans_xyz[1],
																								    ((pUnit->Bond[i].A1Z+pUnit->Bond[i].A2Z)/2) + Trans_xyz[2]);
					    ar->WriteString(szBuffer);
					    sprintf(szBuffer," pigment {color ColorFor%0.2s}}\n",ED[(int)pUnit->Bond[i].A2AtNb-1].Name);
					    ar->WriteString(szBuffer);
				    }else{
					    sprintf(szBuffer,"cylinder{< %f, %f, %f > , < %f, %f, %f>, bondcylRadius",	pUnit->Bond[i].A1X + Trans_xyz[0],
																								    pUnit->Bond[i].A1Y + Trans_xyz[1],
																								    pUnit->Bond[i].A1Z + Trans_xyz[2], 
																								    pUnit->Bond[i].A2X + Trans_xyz[0],
																								    pUnit->Bond[i].A2Y + Trans_xyz[1],
																								    pUnit->Bond[i].A2Z + Trans_xyz[2]);
					    ar->WriteString(szBuffer);
					    ar->WriteString(" pigment {color bondsColor}}\n");
				    }	
                }
			}
		}else{
            std::vector<Molecule> UCMolecules = pUnit->getUnitcell_molecule_cart();
            for(int j=0;j<UCMolecules.size();j++){
                for(int i=0;i< UCMolecules[j].Bond.size();i++){
                    if(UnitCell::Distance(UCMolecules[j].Bond[i].A1X + Trans_xyz[0],
                                          UCMolecules[j].Bond[i].A1Y + Trans_xyz[1],
                                          UCMolecules[j].Bond[i].A1Z + Trans_xyz[2], 
                                          UCMolecules[j].Bond[i].A2X + Trans_xyz[0],
                                          UCMolecules[j].Bond[i].A2Y + Trans_xyz[1],
                                          UCMolecules[j].Bond[i].A2Z + Trans_xyz[2]) > 0.01) {
					    if(DSS.iBondsColorStyle==TwoColors)
					    {
						    sprintf(szBuffer,"cylinder{< %f, %f, %f > , < %f, %f, %f>, ",	UCMolecules[j].Bond[i].A1X + Trans_xyz[0],
																						    UCMolecules[j].Bond[i].A1Y + Trans_xyz[1],
																						    UCMolecules[j].Bond[i].A1Z + Trans_xyz[2], 
																						    ((UCMolecules[j].Bond[i].A1X+UCMolecules[j].Bond[i].A2X)/2) + Trans_xyz[0],
																						    ((UCMolecules[j].Bond[i].A1Y+UCMolecules[j].Bond[i].A2Y)/2) + Trans_xyz[1],
																						    ((UCMolecules[j].Bond[i].A1Z+UCMolecules[j].Bond[i].A2Z)/2) + Trans_xyz[2]);
						    ar->WriteString(szBuffer);
						    ar->WriteString("bondcylRadius");
						    sprintf(szBuffer," pigment {color ColorFor%0.2s}}\n",ED[(int)UCMolecules[j].Bond[i].A1AtNb-1].Name);
						    ar->WriteString(szBuffer);

						    sprintf(szBuffer,"cylinder{< %f, %f, %f > , < %f, %f, %f>, bondcylRadius",	UCMolecules[j].Bond[i].A2X + Trans_xyz[0],
																									    UCMolecules[j].Bond[i].A2Y + Trans_xyz[1],
																									    UCMolecules[j].Bond[i].A2Z + Trans_xyz[2], 
																									    ((UCMolecules[j].Bond[i].A1X+UCMolecules[j].Bond[i].A2X)/2) + Trans_xyz[0],
																									    ((UCMolecules[j].Bond[i].A1Y+UCMolecules[j].Bond[i].A2Y)/2) + Trans_xyz[1],
																									    ((UCMolecules[j].Bond[i].A1Z+UCMolecules[j].Bond[i].A2Z)/2) + Trans_xyz[2]);
						    ar->WriteString(szBuffer);
						    sprintf(szBuffer," pigment {color ColorFor%0.2s}}\n",ED[(int)UCMolecules[j].Bond[i].A2AtNb-1].Name);
						    ar->WriteString(szBuffer);
					    }else{
						    sprintf(szBuffer,"cylinder{< %f, %f, %f > , < %f, %f, %f>, bondcylRadius",	UCMolecules[j].Bond[i].A1X + Trans_xyz[0],
																									    UCMolecules[j].Bond[i].A1Y + Trans_xyz[1],
																									    UCMolecules[j].Bond[i].A1Z + Trans_xyz[2], 
																									    UCMolecules[j].Bond[i].A2X + Trans_xyz[0],
																									    UCMolecules[j].Bond[i].A2Y + Trans_xyz[1],
																									    UCMolecules[j].Bond[i].A2Z + Trans_xyz[2]);
						    ar->WriteString(szBuffer);
						    ar->WriteString(" pigment {color bondsColor}}\n");
					    }	
                    }
				}
			}
		}
 
	//output cell
	//x
	if(DSS.bUnitCellVisible==TRUE){
		sprintf(szBuffer,"cylinder{< %f, %f, %f > , < %f, %f, %f>, cellcylRadius",	pUnit->GetUnitPoint_cart(0).X + Trans_xyz[0],
																					pUnit->GetUnitPoint_cart(0).Y + Trans_xyz[1],
																					pUnit->GetUnitPoint_cart(0).Z + Trans_xyz[2], 
																					pUnit->GetUnitPoint_cart(1).X + Trans_xyz[0],
																					pUnit->GetUnitPoint_cart(1).Y + Trans_xyz[1],
																					pUnit->GetUnitPoint_cart(1).Z + Trans_xyz[2]);
		ar->WriteString(szBuffer);
		ar->WriteString(" pigment {color rgb <1, 0, 0>}}\n");
		//y
		sprintf(szBuffer,"cylinder{< %f, %f, %f > , < %f, %f, %f>, cellcylRadius",	pUnit->GetUnitPoint_cart(0).X + Trans_xyz[0],
																					pUnit->GetUnitPoint_cart(0).Y + Trans_xyz[1],
																					pUnit->GetUnitPoint_cart(0).Z + Trans_xyz[2],
																					pUnit->GetUnitPoint_cart(2).X + Trans_xyz[0],
																					pUnit->GetUnitPoint_cart(2).Y + Trans_xyz[1],
																					pUnit->GetUnitPoint_cart(2).Z + Trans_xyz[2]);
		ar->WriteString(szBuffer);
		ar->WriteString(" pigment {color rgb <0, 1, 0>}}\n");
		sprintf(szBuffer,"cylinder{< %f, %f, %f > , < %f, %f, %f>, cellcylRadius",	pUnit->GetUnitPoint_cart(3).X + Trans_xyz[0],
																					pUnit->GetUnitPoint_cart(3).Y + Trans_xyz[1],
																					pUnit->GetUnitPoint_cart(3).Z + Trans_xyz[2], 
																					pUnit->GetUnitPoint_cart(1).X + Trans_xyz[0],
																					pUnit->GetUnitPoint_cart(1).Y + Trans_xyz[1],
																					pUnit->GetUnitPoint_cart(1).Z + Trans_xyz[2]);
		ar->WriteString(szBuffer);
		ar->WriteString(" pigment {color rgb <0.8, 0.8, 0.8>}}\n");
		sprintf(szBuffer,"cylinder{< %f, %f, %f > , < %f, %f, %f>, cellcylRadius",	pUnit->GetUnitPoint_cart(2).X + Trans_xyz[0],
																					pUnit->GetUnitPoint_cart(2).Y + Trans_xyz[1],
																					pUnit->GetUnitPoint_cart(2).Z + Trans_xyz[2], 
																					pUnit->GetUnitPoint_cart(3).X + Trans_xyz[0],
																					pUnit->GetUnitPoint_cart(3).Y + Trans_xyz[1],
																					pUnit->GetUnitPoint_cart(3).Z + Trans_xyz[2]);
		ar->WriteString(szBuffer);
		ar->WriteString(" pigment {color rgb <0.8, 0.8, 0.8>}}\n");

		sprintf(szBuffer,"cylinder{< %f, %f, %f > , < %f, %f, %f>, cellcylRadius",	pUnit->GetUnitPoint_cart(0).X + Trans_xyz[0],
																					pUnit->GetUnitPoint_cart(0).Y + Trans_xyz[1],
																					pUnit->GetUnitPoint_cart(0).Z + Trans_xyz[2], 
																					pUnit->GetUnitPoint_cart(4).X + Trans_xyz[0],
																					pUnit->GetUnitPoint_cart(4).Y + Trans_xyz[1],
																					pUnit->GetUnitPoint_cart(4).Z + Trans_xyz[2]);
		ar->WriteString(szBuffer);
		ar->WriteString(" pigment {color rgb <0, 0, 1>}}\n");
		sprintf(szBuffer,"cylinder{< %f, %f, %f > , < %f, %f, %f>, cellcylRadius",	pUnit->GetUnitPoint_cart(5).X + Trans_xyz[0],
																					pUnit->GetUnitPoint_cart(5).Y + Trans_xyz[1],
																					pUnit->GetUnitPoint_cart(5).Z + Trans_xyz[2], 
																					pUnit->GetUnitPoint_cart(1).X + Trans_xyz[0],
																					pUnit->GetUnitPoint_cart(1).Y + Trans_xyz[1],
																					pUnit->GetUnitPoint_cart(1).Z + Trans_xyz[2]);
		ar->WriteString(szBuffer);
		ar->WriteString(" pigment {color rgb <0.8, 0.8, 0.8>}}\n");
		sprintf(szBuffer,"cylinder{< %f, %f, %f > , < %f, %f, %f>, cellcylRadius",	pUnit->GetUnitPoint_cart(2).X + Trans_xyz[0],
																					pUnit->GetUnitPoint_cart(2).Y + Trans_xyz[1],
																					pUnit->GetUnitPoint_cart(2).Z + Trans_xyz[2], 
																					pUnit->GetUnitPoint_cart(6).X + Trans_xyz[0],
																					pUnit->GetUnitPoint_cart(6).Y + Trans_xyz[1],
																					pUnit->GetUnitPoint_cart(6).Z + Trans_xyz[2]);
		ar->WriteString(szBuffer);
		ar->WriteString(" pigment {color rgb <0.8, 0.8, 0.8>}}\n");
		sprintf(szBuffer,"cylinder{< %f, %f, %f > , < %f, %f, %f>, cellcylRadius",	pUnit->GetUnitPoint_cart(3).X + Trans_xyz[0],
																					pUnit->GetUnitPoint_cart(3).Y + Trans_xyz[1],
																					pUnit->GetUnitPoint_cart(3).Z + Trans_xyz[2], 
																					pUnit->GetUnitPoint_cart(7).X + Trans_xyz[0],
																					pUnit->GetUnitPoint_cart(7).Y + Trans_xyz[1],
																					pUnit->GetUnitPoint_cart(7).Z + Trans_xyz[2]);
		ar->WriteString(szBuffer);
		ar->WriteString(" pigment {color rgb <0.8, 0.8, 0.8>}}\n");
		sprintf(szBuffer,"cylinder{< %f, %f, %f > , < %f, %f, %f>, cellcylRadius",	pUnit->GetUnitPoint_cart(4).X + Trans_xyz[0],
																					pUnit->GetUnitPoint_cart(4).Y + Trans_xyz[1],
																					pUnit->GetUnitPoint_cart(4).Z + Trans_xyz[2], 
																					pUnit->GetUnitPoint_cart(5).X + Trans_xyz[0],
																					pUnit->GetUnitPoint_cart(5).Y + Trans_xyz[1],
																					pUnit->GetUnitPoint_cart(5).Z + Trans_xyz[2]);
		ar->WriteString(szBuffer);
		ar->WriteString(" pigment {color rgb <0.8, 0.8, 0.8>}}\n");
		sprintf(szBuffer,"cylinder{< %f, %f, %f > , < %f, %f, %f>, cellcylRadius",	pUnit->GetUnitPoint_cart(4).X + Trans_xyz[0],
																					pUnit->GetUnitPoint_cart(4).Y + Trans_xyz[1],
																					pUnit->GetUnitPoint_cart(4).Z + Trans_xyz[2], 
																					pUnit->GetUnitPoint_cart(6).X + Trans_xyz[0],
																					pUnit->GetUnitPoint_cart(6).Y + Trans_xyz[1],
																					pUnit->GetUnitPoint_cart(6).Z + Trans_xyz[2]);
		ar->WriteString(szBuffer);
		ar->WriteString(" pigment {color rgb <0.8, 0.8, 0.8>}}\n");
		sprintf(szBuffer,"cylinder{< %f, %f, %f > , < %f, %f, %f>, cellcylRadius",	pUnit->GetUnitPoint_cart(6).X + Trans_xyz[0],
																					pUnit->GetUnitPoint_cart(6).Y + Trans_xyz[1],
																					pUnit->GetUnitPoint_cart(6).Z + Trans_xyz[2], 
																					pUnit->GetUnitPoint_cart(7).X + Trans_xyz[0],
																					pUnit->GetUnitPoint_cart(7).Y + Trans_xyz[1],
																					pUnit->GetUnitPoint_cart(7).Z + Trans_xyz[2]);
		ar->WriteString(szBuffer);
		ar->WriteString(" pigment {color rgb <0.8, 0.8, 0.8>}}\n");
		sprintf(szBuffer,"cylinder{< %f, %f, %f > , < %f, %f, %f>, cellcylRadius",	pUnit->GetUnitPoint_cart(5).X + Trans_xyz[0],
																					pUnit->GetUnitPoint_cart(5).Y + Trans_xyz[1],
																					pUnit->GetUnitPoint_cart(5).Z + Trans_xyz[2], 
																					pUnit->GetUnitPoint_cart(7).X + Trans_xyz[0],
																					pUnit->GetUnitPoint_cart(7).Y + Trans_xyz[1],
																					pUnit->GetUnitPoint_cart(7).Z + Trans_xyz[2]);
		ar->WriteString(szBuffer);
		ar->WriteString(" pigment {color rgb <0.8, 0.8, 0.8>}}\n");
	}
	
}
void IO::POVDrawMap(VoxelMap *pMap[], int nbMap, CArchive *ar, DisplayStyleStruct DSS, LevelControlSetup LCS)
{//0 Contour lines, 1 Surface, 2 Contour&Surface, 3 None
	float p1[3], p2[3];
	char szBuffer[400];

	for(int i=0;i<nbMap;i++){
		//map as lines
		if(((LCS.iMapMode[i]==0)|(LCS.iMapMode[i]==2))){
			for(int j=0;j<3;j++){
				if(pMap[i]->layer_visibility[j]){
					ar->WriteString("union {");
					for (int k=0;k<pMap[i]->LineMap3D.iNcontourlines[j];k++){ // individual line loop
						p1[0] = pMap[i]->LineMap3D.contourlines[j][k].p[0].x;
						p1[1] = pMap[i]->LineMap3D.contourlines[j][k].p[0].y;
						p1[2] = pMap[i]->LineMap3D.contourlines[j][k].p[0].z;
						p2[0] = pMap[i]->LineMap3D.contourlines[j][k].p[1].x;
						p2[1] = pMap[i]->LineMap3D.contourlines[j][k].p[1].y;
						p2[2] = pMap[i]->LineMap3D.contourlines[j][k].p[1].z;
					
						if(distance(p1,p2)>0.001) // eliminate degenerated cylindres ...
						{
							sprintf(szBuffer,"cylinder{< %f , %f , %f > , ",p1[0],p1[1],p1[2]);
							ar->WriteString(szBuffer);
							sprintf(szBuffer,"< %f , %f , %f > , ",p2[0], p2[1], p2[2]);
							ar->WriteString(szBuffer);
							ar->WriteString("MapcylRadius}\n");
						}
					}
					sprintf(szBuffer,"pigment {color rgb <%f, %f, %f>}}",	pMap[i]->layer_color[j][0],
																		pMap[i]->layer_color[j][1],
																		pMap[i]->layer_color[j][2]);
					ar->WriteString(szBuffer);
				}
			}
		}
		if((LCS.iMapMode[i]==1)|(LCS.iMapMode[i]==2)){
			ar->WriteString("union {");
			for(int j=0;j<pMap[i]->triangles.n;j++){
				ar->WriteString("smooth_triangle{\n");
				sprintf(szBuffer,"<%f, %f, %f>, <%f, %f, %f>\n", pMap[i]->grids.point[pMap[i]->triangles.triangl[j].a].absolut.x,
																 pMap[i]->grids.point[pMap[i]->triangles.triangl[j].a].absolut.y,
																 pMap[i]->grids.point[pMap[i]->triangles.triangl[j].a].absolut.z,
																 pMap[i]->grids.point[pMap[i]->triangles.triangl[j].a].normale.x,
																 pMap[i]->grids.point[pMap[i]->triangles.triangl[j].a].normale.y,
																 pMap[i]->grids.point[pMap[i]->triangles.triangl[j].a].normale.z);
				ar->WriteString(szBuffer);

				sprintf(szBuffer,"<%f, %f, %f>, <%f, %f, %f>\n", pMap[i]->grids.point[pMap[i]->triangles.triangl[j].b].absolut.x,
																 pMap[i]->grids.point[pMap[i]->triangles.triangl[j].b].absolut.y,
																 pMap[i]->grids.point[pMap[i]->triangles.triangl[j].b].absolut.z,
																 pMap[i]->grids.point[pMap[i]->triangles.triangl[j].b].normale.x,
																 pMap[i]->grids.point[pMap[i]->triangles.triangl[j].b].normale.y,
																 pMap[i]->grids.point[pMap[i]->triangles.triangl[j].b].normale.z);
				ar->WriteString(szBuffer);

				sprintf(szBuffer,"<%f, %f, %f>, <%f, %f, %f>", pMap[i]->grids.point[pMap[i]->triangles.triangl[j].c].absolut.x,
																 pMap[i]->grids.point[pMap[i]->triangles.triangl[j].c].absolut.y,
																 pMap[i]->grids.point[pMap[i]->triangles.triangl[j].c].absolut.z,
																 pMap[i]->grids.point[pMap[i]->triangles.triangl[j].c].normale.x,
																 pMap[i]->grids.point[pMap[i]->triangles.triangl[j].c].normale.y,
																 pMap[i]->grids.point[pMap[i]->triangles.triangl[j].c].normale.z);
				ar->WriteString(szBuffer);
				ar->WriteString("}\n");
				
			}
			sprintf(szBuffer,"pigment{color Map%dLevel0Color}}\n",i); 
			ar->WriteString(szBuffer);
		}

	}	
}
void IO::VRMLDrawUnitCell(UnitCell *pUnit, VoxelMap *pMap, int x, int y, int z, CArchive *ar, DisplayStyleStruct DSS, ElementDescription ED[255])
{
	float Trans_xyz[3];
	char szBuffer[400];

	//compute real vector from x, y, z
	Trans_xyz[0] = (float) x * (pUnit->GetUnitPoint_cart(1).X - pUnit->GetUnitPoint_cart(0).X)
				 + (float) y * (pUnit->GetUnitPoint_cart(2).X - pUnit->GetUnitPoint_cart(0).X)
				 + (float) z * (pUnit->GetUnitPoint_cart(4).X - pUnit->GetUnitPoint_cart(0).X);
	Trans_xyz[1] = (float) x * (pUnit->GetUnitPoint_cart(1).Y - pUnit->GetUnitPoint_cart(0).Y) 
				 + (float) y * (pUnit->GetUnitPoint_cart(2).Y - pUnit->GetUnitPoint_cart(0).Y)
				 + (float) z * (pUnit->GetUnitPoint_cart(4).Y - pUnit->GetUnitPoint_cart(0).Y);
	Trans_xyz[2] = (float) x * (pUnit->GetUnitPoint_cart(1).Z - pUnit->GetUnitPoint_cart(0).Z) 
				 + (float) y * (pUnit->GetUnitPoint_cart(2).Z - pUnit->GetUnitPoint_cart(0).Z) 
				 + (float) z * (pUnit->GetUnitPoint_cart(4).Z - pUnit->GetUnitPoint_cart(0).Z);
	
	//output atoms
    std::vector<Atom> unitcell_atoms = pUnit->getUnitcell_atoms_cart();
	if((DSS.DisplayStyle!=Off)&&(DSS.DisplayStyle!=Line)) {
        for(int i=0;i<unitcell_atoms.size();i++) {

            sprintf(szBuffer,"USE mat_%0.2s\n",unitcell_atoms[i].AtomType.MakeUpper());
		    ar->WriteString(szBuffer);
		    sprintf(szBuffer,"Separator {\nTransform {\n");
		    ar->WriteString(szBuffer);
		    sprintf(szBuffer,"translation %f %f %f\n",	unitcell_atoms[i].X + Trans_xyz[0],
													    unitcell_atoms[i].Y + Trans_xyz[1],
													    unitcell_atoms[i].Z + Trans_xyz[2]);
		    ar->WriteString(szBuffer);

		    if ((DSS.DisplayStyle==Ball_and_line)|(DSS.DisplayStyle==Ball_and_cylinder))
			    sprintf(szBuffer,"scaleFactor %f %f %f\n", DSS.BallSize, DSS.BallSize, DSS.BallSize);
		    if(DSS.DisplayStyle==Cylinder)
			    sprintf(szBuffer,"scaleFactor %f %f %f\n", DSS.fCylinderSize, DSS.fCylinderSize, DSS.fCylinderSize);
		    if (DSS.DisplayStyle==Scaled_ball_and_cylinder)
			    sprintf(szBuffer,"scaleFactor %f %f %f\n", DSS.BallSize*DSS.ScaledFactor, DSS.BallSize*DSS.ScaledFactor, DSS.BallSize*DSS.ScaledFactor);

		    ar->WriteString(szBuffer);

		    if (i==0) ar->WriteString("}\nDEF sphere Sphere {}\n}\n\n");
		    else ar->WriteString("}\nUSE sphere\n}\n\n");
        }
	}

	//output bonds
	bool first = true;
	CString tmp;
    if(!DSS.bOrganic) {		
        for(int i=0;i<pUnit->Bond.size();i++) {     
            if(DSS.iBondsColorStyle==TwoColors) {
                sprintf(szBuffer,"mat_%0.2s",ED[(int)pUnit->Bond[i].A1AtNb-1].Name.MakeUpper());
			    tmp = szBuffer;
			    VRMLSaveBond(ar,pUnit->Bond[i].A1X + Trans_xyz[0],
							    pUnit->Bond[i].A1Y + Trans_xyz[1],
							    pUnit->Bond[i].A1Z + Trans_xyz[2], 
							    ((pUnit->Bond[i].A1X+pUnit->Bond[i].A2X)/2) + Trans_xyz[0],
							    ((pUnit->Bond[i].A1Y+pUnit->Bond[i].A2Y)/2) + Trans_xyz[1],
							    ((pUnit->Bond[i].A1Z+pUnit->Bond[i].A2Z)/2) + Trans_xyz[2],
							    DSS.fCylinderSize,
							    first,
							    tmp);
			    first=false;
			    sprintf(szBuffer,"mat_%0.2s",ED[(int)pUnit->Bond[i].A2AtNb-1].Name.MakeUpper());
			    tmp = szBuffer;
			    VRMLSaveBond(ar,pUnit->Bond[i].A2X + Trans_xyz[0],
							    pUnit->Bond[i].A2Y + Trans_xyz[1],
							    pUnit->Bond[i].A2Z + Trans_xyz[2], 
							    ((pUnit->Bond[i].A1X+pUnit->Bond[i].A2X)/2) + Trans_xyz[0],
							    ((pUnit->Bond[i].A1Y+pUnit->Bond[i].A2Y)/2) + Trans_xyz[1],
							    ((pUnit->Bond[i].A1Z+pUnit->Bond[i].A2Z)/2) + Trans_xyz[2],
							    DSS.fCylinderSize,
							    first,
							    tmp);
            } else {
			    tmp = "mat_bond";
			    VRMLSaveBond(ar,pUnit->Bond[i].A1X + Trans_xyz[0],
							pUnit->Bond[i].A1Y + Trans_xyz[1],
							pUnit->Bond[i].A1Z + Trans_xyz[2], 
							pUnit->Bond[i].A2X + Trans_xyz[0],
							pUnit->Bond[i].A2Y + Trans_xyz[1],
							pUnit->Bond[i].A2Z + Trans_xyz[2],
							DSS.fCylinderSize,
							first,
							tmp);
                first=false;
           
            }
        }
	} else {
		Atom *tmpAt1 = new Atom("");
		Atom *tmpAt2 = new Atom("");
		float color1[3], color2[3];
        std::vector<Molecule> UCMolecules = pUnit->getUnitcell_molecule_cart();
        for(int i=0;i<UCMolecules.size();i++) {
            for(int j=0;j<UCMolecules[i].Bond.size();j++)
			{
                 if(DSS.iBondsColorStyle==TwoColors) {
                    sprintf(szBuffer,"mat_%0.2s",ED[(int)UCMolecules[i].Bond[j].A1AtNb-1].Name.MakeUpper());
			        tmp = szBuffer;
			        VRMLSaveBond(ar,UCMolecules[i].Bond[j].A1X + Trans_xyz[0],
							        UCMolecules[i].Bond[j].A1Y + Trans_xyz[1],
							        UCMolecules[i].Bond[j].A1Z + Trans_xyz[2], 
							        ((UCMolecules[i].Bond[j].A1X+UCMolecules[i].Bond[j].A2X)/2) + Trans_xyz[0],
							        ((UCMolecules[i].Bond[j].A1Y+UCMolecules[i].Bond[j].A2Y)/2) + Trans_xyz[1],
							        ((UCMolecules[i].Bond[j].A1Z+UCMolecules[i].Bond[j].A2Z)/2) + Trans_xyz[2],
							        DSS.fCylinderSize,
							        first,
							        tmp);
			        first=false;
			        sprintf(szBuffer,"mat_%0.2s",ED[(int)UCMolecules[i].Bond[j].A2AtNb-1].Name.MakeUpper());
			        tmp = szBuffer;
			        VRMLSaveBond(ar,UCMolecules[i].Bond[j].A2X + Trans_xyz[0],
							        UCMolecules[i].Bond[j].A2Y + Trans_xyz[1],
							        UCMolecules[i].Bond[j].A2Z + Trans_xyz[2], 
							        ((UCMolecules[i].Bond[j].A1X+UCMolecules[i].Bond[j].A2X)/2) + Trans_xyz[0],
							        ((UCMolecules[i].Bond[j].A1Y+UCMolecules[i].Bond[j].A2Y)/2) + Trans_xyz[1],
							        ((UCMolecules[i].Bond[j].A1Z+UCMolecules[i].Bond[j].A2Z)/2) + Trans_xyz[2],
							        DSS.fCylinderSize,
							        first,
							        tmp);
                } else {
			        tmp = "mat_bond";
			        VRMLSaveBond(ar,UCMolecules[i].Bond[j].A1X + Trans_xyz[0],
							    UCMolecules[i].Bond[j].A1Y + Trans_xyz[1],
							    UCMolecules[i].Bond[j].A1Z + Trans_xyz[2], 
							    UCMolecules[i].Bond[j].A2X + Trans_xyz[0],
							    UCMolecules[i].Bond[j].A2Y + Trans_xyz[1],
							    UCMolecules[i].Bond[j].A2Z + Trans_xyz[2],
							    DSS.fCylinderSize,
							    first,
							    tmp);
                    first=false;
           
                }
			}
		}
	}
	//output cell
	VRMLSaveBond(ar,pUnit->GetUnitPoint_cart(0).X + Trans_xyz[0],
					pUnit->GetUnitPoint_cart(0).Y + Trans_xyz[1],
					pUnit->GetUnitPoint_cart(0).Z + Trans_xyz[2], 
					pUnit->GetUnitPoint_cart(1).X + Trans_xyz[0],
					pUnit->GetUnitPoint_cart(1).Y + Trans_xyz[1],
					pUnit->GetUnitPoint_cart(1).Z + Trans_xyz[2],
					DSS.fUnitCellLineSize/100,
					false,
					"mat_Cell");
	//y
	VRMLSaveBond(ar,pUnit->GetUnitPoint_cart(0).X + Trans_xyz[0],
					pUnit->GetUnitPoint_cart(0).Y + Trans_xyz[1],
					pUnit->GetUnitPoint_cart(0).Z + Trans_xyz[2],
					pUnit->GetUnitPoint_cart(2).X + Trans_xyz[0],
					pUnit->GetUnitPoint_cart(2).Y + Trans_xyz[1],
					pUnit->GetUnitPoint_cart(2).Z + Trans_xyz[2], 
					DSS.fUnitCellLineSize/100,
					false,
					"mat_Cell");
	VRMLSaveBond(ar,pUnit->GetUnitPoint_cart(3).X + Trans_xyz[0],
					pUnit->GetUnitPoint_cart(3).Y + Trans_xyz[1],
					pUnit->GetUnitPoint_cart(3).Z + Trans_xyz[2], 
					pUnit->GetUnitPoint_cart(1).X + Trans_xyz[0],
					pUnit->GetUnitPoint_cart(1).Y + Trans_xyz[1],
					pUnit->GetUnitPoint_cart(1).Z + Trans_xyz[2], 
					DSS.fUnitCellLineSize/100,
					false,
					"mat_Cell");
	VRMLSaveBond(ar,pUnit->GetUnitPoint_cart(2).X + Trans_xyz[0],
					pUnit->GetUnitPoint_cart(2).Y + Trans_xyz[1],
					pUnit->GetUnitPoint_cart(2).Z + Trans_xyz[2], 
					pUnit->GetUnitPoint_cart(3).X + Trans_xyz[0],
					pUnit->GetUnitPoint_cart(3).Y + Trans_xyz[1],
					pUnit->GetUnitPoint_cart(3).Z + Trans_xyz[2], 
					DSS.fUnitCellLineSize/100,
					false,
					"mat_Cell");

	VRMLSaveBond(ar,pUnit->GetUnitPoint_cart(0).X + Trans_xyz[0],
					pUnit->GetUnitPoint_cart(0).Y + Trans_xyz[1],
					pUnit->GetUnitPoint_cart(0).Z + Trans_xyz[2], 
					pUnit->GetUnitPoint_cart(4).X + Trans_xyz[0],
					pUnit->GetUnitPoint_cart(4).Y + Trans_xyz[1],
					pUnit->GetUnitPoint_cart(4).Z + Trans_xyz[2], 
					DSS.fUnitCellLineSize/100,
					false,
					"mat_Cell");
	VRMLSaveBond(ar,pUnit->GetUnitPoint_cart(5).X + Trans_xyz[0],
					pUnit->GetUnitPoint_cart(5).Y + Trans_xyz[1],
					pUnit->GetUnitPoint_cart(5).Z + Trans_xyz[2], 
					pUnit->GetUnitPoint_cart(1).X + Trans_xyz[0],
					pUnit->GetUnitPoint_cart(1).Y + Trans_xyz[1],
					pUnit->GetUnitPoint_cart(1).Z + Trans_xyz[2], 
					DSS.fUnitCellLineSize/100,
					false,
					"mat_Cell");
	VRMLSaveBond(ar,pUnit->GetUnitPoint_cart(2).X + Trans_xyz[0],
					pUnit->GetUnitPoint_cart(2).Y + Trans_xyz[1],
					pUnit->GetUnitPoint_cart(2).Z + Trans_xyz[2], 
					pUnit->GetUnitPoint_cart(6).X + Trans_xyz[0],
					pUnit->GetUnitPoint_cart(6).Y + Trans_xyz[1],
					pUnit->GetUnitPoint_cart(6).Z + Trans_xyz[2], 
					DSS.fUnitCellLineSize/100,
					false,
					"mat_Cell");
	VRMLSaveBond(ar,pUnit->GetUnitPoint_cart(3).X + Trans_xyz[0],
					pUnit->GetUnitPoint_cart(3).Y + Trans_xyz[1],
					pUnit->GetUnitPoint_cart(3).Z + Trans_xyz[2], 
					pUnit->GetUnitPoint_cart(7).X + Trans_xyz[0],
					pUnit->GetUnitPoint_cart(7).Y + Trans_xyz[1],
					pUnit->GetUnitPoint_cart(7).Z + Trans_xyz[2], 
					DSS.fUnitCellLineSize/100,
					false,
					"mat_Cell");
	VRMLSaveBond(ar,pUnit->GetUnitPoint_cart(4).X + Trans_xyz[0],
					pUnit->GetUnitPoint_cart(4).Y + Trans_xyz[1],
					pUnit->GetUnitPoint_cart(4).Z + Trans_xyz[2], 
					pUnit->GetUnitPoint_cart(5).X + Trans_xyz[0],
					pUnit->GetUnitPoint_cart(5).Y + Trans_xyz[1],
					pUnit->GetUnitPoint_cart(5).Z + Trans_xyz[2], 
					DSS.fUnitCellLineSize/100,
					false,
					"mat_Cell");
	VRMLSaveBond(ar,pUnit->GetUnitPoint_cart(4).X + Trans_xyz[0],
					pUnit->GetUnitPoint_cart(4).Y + Trans_xyz[1],
					pUnit->GetUnitPoint_cart(4).Z + Trans_xyz[2], 
					pUnit->GetUnitPoint_cart(6).X + Trans_xyz[0],
					pUnit->GetUnitPoint_cart(6).Y + Trans_xyz[1],
					pUnit->GetUnitPoint_cart(6).Z + Trans_xyz[2], 
					DSS.fUnitCellLineSize/100,
					false,
					"mat_Cell");
	VRMLSaveBond(ar,pUnit->GetUnitPoint_cart(6).X + Trans_xyz[0],
					pUnit->GetUnitPoint_cart(6).Y + Trans_xyz[1],
					pUnit->GetUnitPoint_cart(6).Z + Trans_xyz[2], 
					pUnit->GetUnitPoint_cart(7).X + Trans_xyz[0],
					pUnit->GetUnitPoint_cart(7).Y + Trans_xyz[1],
					pUnit->GetUnitPoint_cart(7).Z + Trans_xyz[2], 
					DSS.fUnitCellLineSize/100,
					false,
					"mat_Cell");
	VRMLSaveBond(ar,pUnit->GetUnitPoint_cart(5).X + Trans_xyz[0],
					pUnit->GetUnitPoint_cart(5).Y + Trans_xyz[1],
					pUnit->GetUnitPoint_cart(5).Z + Trans_xyz[2], 
					pUnit->GetUnitPoint_cart(7).X + Trans_xyz[0],
					pUnit->GetUnitPoint_cart(7).Y + Trans_xyz[1],
					pUnit->GetUnitPoint_cart(7).Z + Trans_xyz[2], 
					DSS.fUnitCellLineSize/100,
					false,
					"mat_Cell");
}
void IO::VRMLDrawMap(VoxelMap *pMap[], int nbMap, CArchive *ar, DisplayStyleStruct DSS, LevelControlSetup LCS)
{//0 Contour lines, 1 Surface, 2 Contour&Surface, 3 None
	float p1[3], p2[3];
	char szBuffer[400];
	CString tmp;

	for(int i=0;i<nbMap;i++) {
		//map as lines
		if(((LCS.iMapMode[i]==0)|(LCS.iMapMode[i]==2))){
			for(int j=0;j<3;j++){
				if(pMap[i]->layer_visibility[j]){
					for (int k=0;k<pMap[i]->LineMap3D.iNcontourlines[j];k++){ // individual line loop
						p1[0] = pMap[i]->LineMap3D.contourlines[j][k].p[0].x;
						p1[1] = pMap[i]->LineMap3D.contourlines[j][k].p[0].y;
						p1[2] = pMap[i]->LineMap3D.contourlines[j][k].p[0].z;
						p2[0] = pMap[i]->LineMap3D.contourlines[j][k].p[1].x;
						p2[1] = pMap[i]->LineMap3D.contourlines[j][k].p[1].y;
						p2[2] = pMap[i]->LineMap3D.contourlines[j][k].p[1].z;
					
						if(distance(p1,p2)>0.001) // eliminate degenerated cylindres ...
						{
							sprintf(szBuffer,"mat_%d%d",i, j);
							tmp = szBuffer;
							
							VRMLSaveBond(ar,p1[0],
											p1[1],
											p1[2], 
											p2[0],
											p2[1],
											p2[2], 
											DSS.fMapLineSize/100,
											false,
											tmp);	
						}
					}
				}
			}
		}
        int reverse;
		if((LCS.iMapMode[i]==1)|(LCS.iMapMode[i]==2)) {          
			if (false) reverse = -1; else reverse = 1;
			//if (j==1) break;//At this time for the first map only
			if (i==0) ar->WriteString("USE mat_00\n");
			if (i==1) ar->WriteString("USE mat_10\n");
			ar->WriteString("# Triangle Mesh definition \n");
			ar->WriteString("Separator {\nCoordinate3 {\npoint [\n");

            for(int j=0;j<pMap[i]->grids.n;j++) {
				sprintf(szBuffer,"	%f %f %f,\n",	pMap[i]->grids.point[j].absolut.x,
													pMap[i]->grids.point[j].absolut.y, 
													pMap[i]->grids.point[j].absolut.z);

				ar->WriteString(szBuffer);
						
			}
			ar->WriteString("]\n}\n\n");
						
			ar->WriteString("# Normal definition \n");
			ar->WriteString("Normal { \n	vector [\n");
			for(int j=0;j<pMap[i]->grids.n;j++) {
				sprintf(szBuffer,"	%f %f %f,\n",	pMap[i]->grids.point[j].normale.x * reverse,
													pMap[i]->grids.point[j].normale.y * reverse, 
													pMap[i]->grids.point[j].normale.z * reverse);

				ar->WriteString(szBuffer);
							
			}
			ar->WriteString("]\n}\n\n");
						
			ar->WriteString("NormalBinding {value PER_VERTEX_INDEXED}\nIndexedFaceSet  {\ncoordIndex [\n");
			for(int j=0;j<pMap[i]->triangles.n;j++) {
				sprintf(szBuffer,"	%d, %d, %d",pMap[i]->triangles.triangl[j].a, pMap[i]->triangles.triangl[j].b, pMap[i]->triangles.triangl[j].c);
				ar->WriteString(szBuffer);
				if(i<pMap[i]->triangles.n-1) ar->WriteString(",-1,\n");
					else ar->WriteString("\n");
			}

			ar->WriteString("]\n}\n}\n\n");

		}
    }	
}
void IO::VRMLSaveBond(CArchive *ar, float x1, float y1, float z1, float x2, float y2, float z2, float radius, BOOL first, CString color)
{
	char szBuffer[400];
	float Sx, Sy, Sz, h;
	float x, y, z;
	
	sprintf(szBuffer,"USE %s\n",color);
	ar->WriteString(szBuffer);
	
	sprintf(szBuffer,"Separator {\nTransform {\n");
	ar->WriteString(szBuffer);

	Sx=(x1+x2)/2;
	Sy=(y1+y2)/2;
	Sz=(z1+z2)/2;
	sprintf(szBuffer,"translation %f %f %f\n", Sx, Sy, Sz);
	ar->WriteString(szBuffer);
	
	x=x2-x1; y=y2-y1; z=z2-z1; // v = (x,y,z) = AB = B-A
	h=(float) sqrt(pow(x,2)+pow(y,2)+pow(z,2)); // |AB|
    if(h!=0) {
	    x/=h;
	    y/=h;
	    z/=h;
    }
	
	z=z+1;

	sprintf(szBuffer,"rotation %f %f %f 3.1413\n}\n", x, y, z);
	ar->WriteString(szBuffer);

	sprintf(szBuffer,"Transform {\n");
	ar->WriteString(szBuffer);

	sprintf(szBuffer,"rotation 1 0 0 1.57\n");
	ar->WriteString(szBuffer);
	
	sprintf(szBuffer,"scaleFactor %f %f %f\n", radius, h, radius);
	ar->WriteString(szBuffer);

	if (first) ar->WriteString("}\nDEF cylinder Cylinder {\nparts SIDES\nheight 1\n}\n}\n\n");
		else ar->WriteString("}\nUSE cylinder\n}\n\n");

}
float IO::distance(float p1[3], float p2[3])
{
// simple distance of two points
float d1,d2,d3;

	d1=p1[0]-p2[0];
	d2=p1[1]-p2[1];
	d3=p1[2]-p2[2];

return (float)sqrt(d1*d1+d2*d2+d3*d3);
}
void IO::writeLog(CString log) 
{
	std::ofstream outfile;
    log += "\n"; 
    outfile.open("log_MCE.txt", std::ios::out | std::ios::app);
	outfile.write(log.GetString(), log.GetLength());
    outfile.close();
}