#include "../interface/environment.h"
#include "../delphi/delphi_constants.h"
#include "../misc/misc_grid.h"
#include "space_templates.h"

#include "space.h"

using namespace std;

void CDelphiSpace::msrf()
{
    delphi_integer itot, vtot;
    string fixyn;
    string vdat1, line[6], fname;
    bool fix;
    SGrid <delphi_real> v1,v2,vxyz;
    SGrid <delphi_integer> iv123;

    delphi_integer mxvtx,i,ia1,ia2,ia3,iate,it,imxvtx;
    delphi_integer ib,iv1,iv2,iv3,j,k,mxtri,ntot2;
    delphi_real aa,area,cc,bb,rad,rad2,ss;
    delphi_real tar,tne4,vmg;

    iv1=0;
    iv2=0;
    iv3=0;
    vtot=0;
    itot=0;
    ntot2=0;

    mxvtx=iBoundNum*2.2;
    mxtri=2*mxvtx;

    //epsdim=iNatom+iNObject+2; replaced in space module semi global

    for(k=1; k<=iGrid.nZ; k++)
    {
        for(j=1; j<=iGrid.nY; j++)
        {
            for(i=1; i<=iGrid.nX; i++)
            {
                //changed from div mod, it should only serve the means
                egrid[i][j][k]=iepsmp[i][j][k]/epsdim;
            }
        }
    }

    for(i=1; i<=iGrid.nX; i++)
    {
        for(j=1; j<=iGrid.nY; j++)
        {
            egrid[i][1][j].nX = egrid[i][1][j].nY;
            egrid[i][j][1].nX = egrid[i][j][1].nZ;
        }
    }

    for(k=2; k<=iGrid.nZ; k++)
    {
        for(j=2; j<=iGrid.nY; j++)
        {
            for(i=2; i<=iGrid.nX; i++)
            {
                iate = 0;
                if (egrid[i][j][k].nX > 0) iate = iate + 1;
                if (egrid[i][j][k].nY > 0) iate = iate + 1;
                if (egrid[i][j][k].nZ > 0) iate = iate + 1;
                if (egrid[i-1][j][k].nX > 0) iate = iate + 1;
                if (egrid[i][j-1][k].nY > 0) iate = iate + 1;
                if (egrid[i][j][k-1].nZ > 0) iate = iate + 1;

                if (iate<=3)
                    egrid[i][j][k].nX = 0;
                else
                    egrid[i][j][k].nX = 1;
            }
        }
    }

    vindx.assign(mxtri+1, sgrid_temp_int);
    vert.assign(mxvtx+1, sgrid_temp_real);

    #ifdef VERBOSE
    cout << "mxtri,mxvtx= " << " " << mxtri << " " << mxvtx << endl;
    #endif

    vdat1="./";

    if (vtot>mxvtx)
    {
        cout <<"vtot = " << vtot << " > mxvtx = " << mxvtx << endl;
        cout <<"increase mxvtx in msrf.f" << endl;
        exit(0);
    }

    for(ib=1; ib<=vtot; ib++)
    {
        vert[ib]=vert[ib]/2.;
    }

    itot = itot/3;

    //fScale boundary grid pointeger positions relative to acc data
    #ifdef VERBOSE
    cout <<"scaling vertices" << endl;
    #endif

    vnorm.assign(vtot+1, sgrid_temp_real);
    vnorm2.assign(vtot+1, sgrid_temp_real);

    //fix holes and make vertex to triangle arrays allocate hole-fixing arrays next hole-fixing variables
    if (vtot < mxvtx/2)
        imxvtx = vtot*2;
    else
        imxvtx = mxvtx;

    vtlen.assign(imxvtx+1,0);
    vtpnt.assign(imxvtx+1,0);
    vtlst.assign(6*imxvtx+1,0);

    get_pt2d <delphi_integer> (tmlst,9+1,imxvtx+1);

    if (ntot2 > 0)
        fix = true;
    else
        fix = false;

    while(fix)
    {
        if (ntot2 > 0)
            fix = true;
        else
            fix = false;
    }

    if (itot>mxtri)
    {
        cout <<"itot = " << itot << " > mxtri = " << mxtri << endl;
        cout <<"increase mxtri in msrf.f" << endl;
        exit(0);
    }

    if(vtlen.size()>0) vector <delphi_integer>().swap(vtlen);
    if(vtlst.size()>0) vector <delphi_integer>().swap(vtlst);
    if(vtpnt.size()>0) vector <delphi_integer>().swap(vtpnt);
    if(tmlst != NULL) free_pt2d(tmlst,9+1,imxvtx+1);


    #ifdef VERBOSE
    cout << "number of vertices = " << vtot << endl;
    cout << "number of triangles = " << itot << endl;
    #endif

    //calculate area
    area=0.0;

    for(it=1; it<=itot; it++)
    {
        iv123=vindx[it];
        v1=vert[iv2]-vert[iv1];
        v2=vert[iv3]-vert[iv1];

        //2011-05-20 Vector product defined in operators_on_coordinates module
        vxyz=optCross(v1,v2);
        vmg=sqrt(optDot(vxyz,vxyz));
        tar=vmg/2.;
        vxyz=vnorm[iv1]+vnorm[iv2]+vnorm[iv3];
        vmg=sqrt(optDot(vxyz,vxyz));

        vnorm2[iv1]=vnorm2[iv1]+(vxyz/vmg);
        vnorm2[iv2]=vnorm2[iv2]+(vxyz/vmg);
        vnorm2[iv3]=vnorm2[iv3]+(vxyz/vmg);

        //calculate spherical triangle area if appropriate
        ia1=atndx[iv1];
        ia2=atndx[iv2];
        ia3=atndx[iv3];

        if (ia1>0)
        {
            if (ia1==ia2&&ia1==ia3)
            {
                rad=sDelPhiPDB[ia1].radius;
                rad2=rad*rad;
                aa=optSum((vert[iv2]-vert[iv1])*(vert[iv2]-vert[iv1]));
                bb=optSum((vert[iv3]-vert[iv2])*(vert[iv3]-vert[iv2]));
                cc=optSum((vert[iv1]-vert[iv3])*(vert[iv1]-vert[iv3]));

                aa=acos(1.-aa/(2.*rad2));
                bb=acos(1.-bb/(2.*rad2));
                cc=acos(1.-cc/(2.*rad2));
                ss=(aa+bb+cc)*.5;
                tne4=sqrt(tan(ss*.5)*tan((ss-aa)*.5)*tan((ss-bb)*.5)*tan((ss-cc)*.5));
                tar=4.*atan(tne4)*rad2;
            }
        }

        area=area+tar;
    }

    for(i=1; i<=vtot; i++)
    {
        vmg=sqrt( optDot( vnorm2[i],vnorm2[i] ) );
        vnorm2[i]=vnorm2[i]/vmg;
    }

    #ifdef VERBOSE
    cout <<"MS area = " << area << endl;
    #endif

    //2011-05-26 Subroutine wrtsurf is short and called only once,
    // therefore put the body of the void here

    //-------- Begin of wrtsurf body -----------------------------

    if (!ibem)
    {
        ofstream surfile;
        surfile.open ("grasp.srf");

        fname="grasp.srf";

        #ifdef VERBOSE
        cout << "writing GRASP file to " << fname << endl;
        #endif

        line[1]="format=2";
        line[2]="vertices,normals,triangles";

        surfile << "format=2" << endl;
        surfile << "vertices,normals,triangles" << endl;
        surfile << endl;
        surfile <<  setw(6) << vtot << setw(6) << itot << setw(6) << iGrid << setw(12) << setprecision(6) << fScale;
        surfile << setw(6) << setw(12) << setprecision(6) << cOldMid.nX
                << setw(6) << setw(12) << setprecision(6) << cOldMid.nY
                << setw(6) << setw(12) << setprecision(6) << cOldMid.nZ;

        #ifdef VERBOSE
        cout << "writing data for" << vtot << " vertices and" << itot << " triangles" << endl;
        #endif

        for(i=1;i<=mxvtx;i++)
            surfile << vert[i] << endl;

        for(i=1;i<=vtot;i++)
            surfile << vnorm[i] << endl;

        for(i=1;i<=mxtri;i++)
            surfile << vindx[i] << endl;

        surfile.close();

        #ifdef VERBOSE
        cout << "finished writing " << fname << endl;
        #endif
    }
    else
    {
        ofstream surfile;
        surfile.open ("bem.srf");

        surfile << vtot << " " << itot << endl;

        for(i=1; i<=vtot; i++)
            surfile << vert[i] << endl;

        for(i=1; i<=itot; i++)
            surfile << vindx[i] << endl;

        for(i=1; i<=vtot; i++)
            surfile << vnorm[i] << endl;

        surfile.close();
    }
    //-------- End of wrtsurf body --------------------------------

}// void msrf;
