// convert ramachandran angle data into a 2d-histogram and
// write it out as a cube file.
// Copyright (c) 2003-2004 by Axel.Kohlmeyer@theochem.ruhr-uni-bochum.de
// $Id: rama2cube.cc,v 1.1 2004/05/16 17:54:14 akohlmey Exp $

#include <iostream>
#include <fstream>
#include <string>
#include <cstdlib>
#include <cstring>
#include <cstdio>

using std::cout;
using std::cin;
using std::istream;
using std::ostream;
using std::ifstream;
using std::ofstream;

int xbins=36;
int ybins=xbins;
int zbins=xbins;

/* strip filename from the full path (similar to basename(1))
 * we accept '/' or '\' as a path separator. a bit ugly, but anybody
 * who uses '/' or '\' in a regular filename deserves to be shot. */
static const char *basename(const char *path)
{
        const char *ptr;

        ptr = strrchr(path, '/');
        if (! ptr) {            /* no '/' path separator */
                ptr = strrchr(path, '\\');
        }

        if (! ptr) {            /* no '\' path separator */
                ptr = path;
        } else {
                ++ ptr; /* skip over separator */
        }
        return ptr;

}


// help message
static void usage(const char *n) 
{
    cerr << "\nusage: \n" << n 
         << " [-x <xbins>] [-y <ybins>] [-z <zbins> ] [-g <descr>] [-i <infile>] [ -o <outfile>]\n";
    cerr << "\n" << n << " converts ramachandran data "
         << "to a gaussian cube file.\n\nOptions:\n";
    cerr << "\t-x\tset number of bins in x-direction (default " << xbins << ")\n";
    cerr << "\t-y\tset number of bins in y-direction (default " << ybins << ")\n";
    cerr << "\t-z\tset number of bins in z-direction (default " << zbins << ")\n";
    cerr << "\t-i\tset input file to <infile>(default stdin)\n";
    cerr << "\t-o\tset output file to <outfile>(default stdout)\n";
}


int main(int argc, char **argv) 
{
    int i,j,k;
    istream *in;
    ostream *out;

    in = &cin;
    out = &cout;
    
// parse commandline
    i = 1;
    while (i < argc) {

        if (string("-h") == argv[i]) {
            usage(basename(argv[0]));
            return 1;
        }


        if (string("-x") == argv[i]) {
            ++i;
            if (i < argc) {
                xbins = atoi(argv[i]);
                if (xbins < 10) {
                    cerr << "xbins (=" << xbins << ") smaller than 10 makes no sense\n";
                    cerr << "try again\n";
                    return 10;
                }
            } else {
                cerr << "no argument to option '-x' \n";
                return 1;
            }
            ++ i;
            continue;
        }

        if (string("-y") == argv[i]) {
            ++i;
            if (i < argc) {
                ybins = atoi(argv[i]);
                if (ybins < 10) {
                    cerr << "ybins (=" << ybins << ") smaller than 10 makes no sense\n";
                    cerr << "try again\n";
                    return 20;
                }
            } else {
                cerr << "no argument to option '-y' \n";
                return 1;
            }
            ++ i;
            continue;
        }

        if (string("-z") == argv[i]) {
            ++i;
            if (i < argc) {
                zbins = atoi(argv[i]);
                if (zbins < 10) {
                    cerr << "zbins (=" << zbins << ") smaller than 10 makes no sense\n";
                    cerr << "try again\n";
                    return 20;
                }
            } else {
                cerr << "no argument to option '-z' \n";
                return 1;
            }
            ++ i;
            continue;
        }

        if (string("-i") == argv[i]) {
            ++i;
            if (i < argc) {
                in = new ifstream(argv[i]);
                if (! *in) {
                    cerr << "could not open input file " << argv[i] << " for reading\n";
                    return 21;
                }
            } else {
                cerr << "no argument to option '-i' \n";
                return 1;
            }
            ++ i;
            continue;
        }

        if (string("-o") == argv[i]) {
            ++i;
            if (i < argc) {
                out = new ofstream(argv[i]);
                if (! *out) {
                    cerr << "could not open output file " << argv[i] << " for writing\n";
                    return 22;
                }
            } else {
                cerr << "no argument to option '-o' \n";
                return 1;
            }
            ++ i;
            continue;
        }

        cerr << "unknown argument: " << argv[i] << endl;
        usage(basename(argv[0]));
        return 2;
    }

    // prepare histogram
    
    int hist[xbins][ybins];
    double x, y, dx(360.0/double(xbins-1)), dy(360.0/double(ybins-1));
    
    for (i = 0; i < xbins; ++i) {
        for (j = 0; j < xbins; ++j) {
            hist[i][j]= 0;
        }
    }

    while (*in) {
        char buf[512];
        in->getline(buf, 512, '\n');
        
        // skip over comments and plot commands
        if (buf[0] == '@') continue;
        if (buf[0] == '#') continue;

        // for genrama.tcl output
        sscanf(buf, "%*d %*s %lg %lg", &x, &y);

        // for gromacs output.
        // sscanf(buf, "%lg %lg", &x, &y);

        i = int((x + 180.0) / dx);
        j = int((y + 180.0) / dy);
        if ((i < 0) || (j < 0) || (i >= xbins) || (j >= ybins)) {
            cerr << "Overflow: x=" << x << "  y=" << y 
                 << " => i=" << i << " j=" << j << endl;
            continue;
        }
        
        ++ hist[i][j];
    }

    // determine scale.
    int zmax = 0;
    // fix histogram
    for (i = 0; i < xbins; ++i) {
        hist[i][ybins-1] = hist[i][0];
    }
    for (j = 0; j < ybins; ++j) {
        hist[xbins-1][j] = hist[0][j];
    }

    for (i = 0; i < xbins; ++i) {
        for (j = 0; j < ybins; ++j) {
            zmax = (zmax < hist[i][j]) ?  hist[i][j] : zmax;
        }
    }

    // now the formatted output madness starts.
    // i really should have written this in fortran...
    *out << " RAMACHANDRAN data cube\n";
    *out << " Total SCF Density\n";
    char outbuf[255];
    sprintf(outbuf, "%5d %12.6f %12.6f %12.6f\n",     8, 0.0, 0.0, 0.0);
    *out << outbuf;
    sprintf(outbuf, "%5d %12.6f %12.6f %12.6f\n", xbins, 0.1, 0.0, 0.0);
    *out << outbuf;
    sprintf(outbuf, "%5d %12.6f %12.6f %12.6f\n", ybins, 0.0, 0.1, 0.0);
    *out << outbuf;
    sprintf(outbuf, "%5d %12.6f %12.6f %12.6f\n", zbins, 0.0, 0.0, 0.03333);
    *out << outbuf;

    sprintf(outbuf, "%5d %12.6f %12.6f %12.6f %12.6f\n", 1, 1.0, double(xbins)*0.05, 0.0, 0.0);
    *out << outbuf;
    sprintf(outbuf, "%5d %12.6f %12.6f %12.6f %12.6f\n", 1, 1.0, 0.0, double(ybins)*0.05, 0.0);
    *out << outbuf;
    sprintf(outbuf, "%5d %12.6f %12.6f %12.6f %12.6f\n", 1, 1.0, double(xbins)*0.05, double(ybins)*0.1, 0.0);
    *out << outbuf;
    sprintf(outbuf, "%5d %12.6f %12.6f %12.6f %12.6f\n", 1, 1.0, double(xbins)*0.1, double(ybins)*0.05, 0.0);
    *out << outbuf;

    sprintf(outbuf, "%5d %12.6f %12.6f %12.6f %12.6f\n", 3, 3.0, 0.0, 0.0, 0.0);
    *out << outbuf;
    sprintf(outbuf, "%5d %12.6f %12.6f %12.6f %12.6f\n", 3, 3.0, 0.0, double(ybins)*0.1, 0.0);
    *out << outbuf;
    sprintf(outbuf, "%5d %12.6f %12.6f %12.6f %12.6f\n", 3, 3.0, double(xbins)*0.1, 0.0, 0.0);
    *out << outbuf;
    sprintf(outbuf, "%5d %12.6f %12.6f %12.6f %12.6f\n", 3, 3.0, double(xbins)*0.1, double(ybins)*0.1, 0.0);
    *out << outbuf;

    
    // now write out the cube data, fortran style: 6 columns maximum,
    // start with a new line for each block of z-values.
    int cube[xbins][ybins][zbins];
    double dz(double(zmax)/double(zbins-1));
    for (i = 0; i < xbins; ++i) {
        for (j = 0; j < ybins; ++j) {
            for (k = 0; k < zbins; ++k) {
                cube[i][j][k] = 0;
            }
            cube[i][j][0] = 1;

            const int zidx = int( double(hist[i][j]) / dz );
            for (k = 0; k < zidx; ++k) {
                cube[i][j][k] = zidx-k;
            }

            for (k = 0; k < zbins; ++k) {
                sprintf(outbuf," %12.5e", 0.01*double(cube[i][j][k]));
                *out << outbuf;
                if ((k+1) % 6 == 0) *out << endl;
            }
            if (k % 6) *out << endl;
        }
    }
    *out << endl;
    
    return 0;
}

    
// Local Variables:
// mode: c++
// compile-command: "g++ -O -Wall -o rama2cube rama2cube.cc"
// End: