#!/usr/bin/perl -w
#
# Script to convert range of z values (center of bilayer becomes 0)
#

use strict;
use warnings;

use List::Util qw( min max );


# Open input files in read mode
open INPUTFILE, "<", $ARGV[0] or die $!;
open INPUTFILE2, "<", $ARGV[1] or die $!;

# Open output file in write mode
open OUTPUTFILE, ">", $ARGV[2] or die $!;

# Read the input file line by line and split in two arrays:
my @z1;
my @density1;
my @z2;
my @density2;

while(my $line = <INPUTFILE>){
    chomp $line;
    my @linearray = split(" ", $line);
#	print "@linearray \t";
    push(@z1, $linearray[0]);
    push(@density1, $linearray[1]);
}
my $numelements1=@z1;
#print "@z1 ";
#print "@density1 ";


my $i = 0;
while(my $line = <INPUTFILE2>){
    chomp $line;
    my @linearray = split(" ", $line);
        #print $line;
	push(@z2, $linearray[0]);
        push(@density2, $linearray[1]);
        $i++;
}
my $numelements2=@density2;

$i = 0;
while($z1[$i] != $z2[0]){
        unshift(@density2,0);
	$i++;
}
while($z1[$i] != $z2[$i]){
        push(@density2,0);
        $numelements2=@density2;
                $i++;
}

my $min = min @z1;
my $max = max @z1;
my $range = $max - $min;
my $diff = $range / 2; 

print OUTPUTFILE "\@type xy \n";
print OUTPUTFILE "\@title \"Density\" \n";
print OUTPUTFILE "\@xaxis label \"z (angstroms)\" \n";
print OUTPUTFILE "\@yaxis label \"Electron density\" \n";
print OUTPUTFILE "\@ s0 legend \"Popc\" \n";
print OUTPUTFILE "\@ s1 legend \"Asp\" \n";
print OUTPUTFILE "\@ legend 0.9, 0.8 \n";
print OUTPUTFILE "\@ world -45 0 45 0.9 \n";
print OUTPUTFILE "\@ WORLD XMIN -45 \n";
print OUTPUTFILE "\@ WORLD XMAX 45 \n";
print OUTPUTFILE "\@ WORLD YMIN 0 \n";
print OUTPUTFILE "\@ WORLD YMAX 0.9 \n";
print OUTPUTFILE "\@ AUTOSCALE ONREAD NONE \n";



for (my $j = 0; $j < $numelements1; $j++){
	$z1[$j] = $j - $diff;
	print OUTPUTFILE "$z1[$j] \t $density1[$j] \t $density2[$j] \n";
}

# Close the input and output files
close INPUTFILE;
close OUTPUTFILE;