#!/usr/bin/perl -w
use Getopt::Std;                # for easy command-line options processing
use File::Copy;                 # File copy

# Command line option processing
getopts("cobdrpas:i:x:y:z:");
$opt_o ||= 0;                   # original coordinates
$opt_d ||= 0;                   # if to display the image
$opt_b ||= 0;                   # display ball-and-stick model
$opt_r ||= 0;                   # only P atoms + base-Ring atoms
$opt_p ||= 0;                   # Protein atoms
$opt_a ||= 0;                   # All atoms
$opt_c ||= 0;                   # if just to clean up temporary files
$opt_s ||= 0;                   # set scale factor: -s=NUM or -s NUM
$opt_s =~ s/=//;
$opt_s = 0 if ($opt_s < 0);
$opt_i ||= 0;                   # output image name
$opt_i =~ s/=//;

$opt_x ||= 0;                   # Rotation axis and angle
$opt_x =~ s/=//;
$opt_y ||= 0;
$opt_y =~ s/=//;
$opt_z ||= 0;
$opt_z =~ s/=//;

# DIRECTORY CONTAINING 3DNA BINARIES: change this line as appropriate
$BDIR = $ENV{X3DNA};
if ($BDIR) {                # X3DNA
    $BDIR =~ s/.$// if (substr($BDIR, -1) eq '/');
    $X3DNA_BIN = "$BDIR/bin";
} else {
    $X3DNA_BIN = "$ENV{HOME}/X3DNA/bin"; # change this line as necessary
}

if (!$opt_c && @ARGV == 1) {
    $pdbfile = $ARGV[0];
    -e $pdbfile or die "PDB file <$pdbfile> does NOT exist!\n";
    if ($pdbfile eq "t.pdb") { # Special case: make a copy
        copy($pdbfile, "__copy__");
        $pdbfile = "__copy__";
    }
}

system("$X3DNA_BIN/dcmnfile -x");
die "\n" if ($opt_c);  # just file clean up

if (@ARGV != 1) {
    print <<BLOCVIEW_USAGE;
===========================================================================
SYNOPSIS
    blocview [OPTION]... PDBFILE
DESCRIPTION
    Generates a schematic image which combines base block representation
    with protein ribbon. The image has informative color coding for the
    nucleic acid part and is set in the "best-view" by default. Users need
    to have MolScript, Raster3D and ImageMagick properly installed on their
    system.
        -o   use original coordinates contained in the PDB data file
        -d   display the generated image using "display" of ImageMagick
        -b   ball and stick model with filled base ring
        -c   clean up temporary common files
        -r   only backbone P atoms + base Ring atoms of nucleic acids
        -p   set the best view based on Protein atoms
        -a   set the best view based on All atoms
        -s[=]NUM    set scale factor for the image
        -i[=]IMAGE  set image file name (default to t.jpg)
        -x|y|z=ANGLE  rotation around x, y, or z-axis by ANGLE degrees
        PDBFILE   a PDB data file name (other than "t.pdb")
EXAMPLES
    blocview -d -i=bdl084.jpg bdl084.pdb
===========================================================================
BLOCVIEW_USAGE
    die "\n";
}

copy($pdbfile, "t.pdb");        # make a copy to a default file name

# Get only protein, NA, ligands and delete H-atoms
system("$X3DNA_BIN/get_part -pnt t.pdb temp");
system("$X3DNA_BIN/get_part -d temp t.pdb");

# Get the best view according to principle moment of inertia
if (!$opt_o) {
    $opt_a = 0 if ($opt_p && $opt_a); # protein preferred than all atoms
    if ($opt_p) {               # protein part
        system("$X3DNA_BIN/rotate_mol -cp t.pdb temp");
    } elsif ($opt_a) {          # all atoms
        system("$X3DNA_BIN/rotate_mol -ca t.pdb temp");
    } else {                    # nucleic acid part
        system("$X3DNA_BIN/rotate_mol -cb t.pdb temp");
    }

    $rotfile = "rotxyz.ang";
    open(FOUT, ">$rotfile") || die "Can't open file <$rotfile> for writing: $!\n";
    print FOUT "by rotation x $opt_x\n";
    print FOUT "by rotation y $opt_y\n";
    print FOUT "by rotation z $opt_z\n";
    close(FOUT);
    system("$X3DNA_BIN/rotate_mol -c -r=$rotfile temp t.pdb");
}

# Get header for the whole structure
if ($opt_s) {                   # user supplied scale: e.g. for biological unit
    system("$X3DNA_BIN/r3d_atom -s=$opt_s t.pdb t.r3d");
} else {                        # whole structure for default scale
    system("$X3DNA_BIN/r3d_atom t.pdb t.r3d");
}
open(FINP, "t.r3d") || die "Can't open file <t.r3d> for reading: $!\n";
open(FOUT, ">t0.r3d") || die "Can't open file <t0.r3d> for writing: $!\n";
while (<FINP>) {
    print FOUT if ($. <= 20);
}
close(FINP);
close(FOUT);

if ($opt_b) {
    copy("t.pdb", "temp");
} else {
    if ($opt_r) {                   # P + base Ring atoms
        system("$X3DNA_BIN/get_part -z t.pdb tb.pdb");
    } else { # Get the nucleic acid backbone (without O1P/O2P) and base ring atoms
        system("$X3DNA_BIN/get_part -x t.pdb tb.pdb");
    }

    # Nucleic acid base and sugar are colored by residue type
    system("$X3DNA_BIN/r3d_atom -ncz -r=0.16 tb.pdb t1.r3d");

    # Block representation of the bases
    system("$X3DNA_BIN/pdb2img -rcnm t.pdb t2.r3d");

    # Use "molauto" to get protein secondary structure & nucleic acid backbone
    system("molauto -nocentre -ss_hb t.pdb > temp"); # H-bond based 2nd str.
    system("$X3DNA_BIN/get_part -c t.pdb temp2");

    # Use "molscript" to get an input to Raster3D, delete its header section
    system("molscript -r < temp2 | del_ms -n > t3.r3d");

    # Get input for ligands
    system("$X3DNA_BIN/get_part -t t.pdb temp");
}

system("$X3DNA_BIN/r3d_atom -r=0.0 -b=0.20 -o -n temp t4.r3d"); # ball for atoms
system("$X3DNA_BIN/r3d_atom -r=0.12 -gn temp t5.r3d"); # gray rod

# Combine all 5 parts together and render it with Raster3D
open(FOUT, ">t.r3d") || die "Cann't open file <t.r3d> for writing: $!\n";
for ($i = 0; $i <= 5; $i++) {
    next if ($opt_b && $i >= 1 & $i <= 3);
    $nfile = "t$i.r3d";
    if (-e $nfile) {
        open(FINP, $nfile) || die "Cann't open file <$nfile> for reading: $!\n";
        print FOUT while (<FINP>);
        close(FINP);
    }
}
system("render < t.r3d > t.avs");

# Crop the image to its minimum size
system("convert -crop 0x0 t.avs t.jpg");

# Save image to its designed name
copy("t.jpg", "$opt_i") if ($opt_i);

# Get special file name back
copy($pdbfile, "t.pdb") if ($pdbfile eq '__copy__');

# Print out scale factor used
open(FINP, "t.r3d") || die "Can't open file <t.r3d> for reading: $!\n";
for ($i = 1; $i <= 16; $i++) {
    $str = <FINP>;
}
close(FINP);

@four = split(" ", $str);
printf "***** Scale factor used: ===> %.2f <===\n", $four[3];

# Finally display the image if -d option is given
system("display t.jpg") if ($opt_d);