#!/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 INPUTFILE3, "<", $ARGV[2] or die $!;
open INPUTFILE4, "<", $ARGV[3] or die $!;

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

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

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

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

my $i = 0;
while(my $line = <INPUTFILE3>){
    chomp $line;
    my @linearray = split(" ", $line);
	push(@z3, $linearray[0]);
	push(@density3, $linearray[1]);
	if ($z1[$i] != $z3[$i]){
		unshift(@z3,$z1[$i]); #add element to TOP of array
		unshift(@density3,0);
	}
	$i++;
}
my $numelements3=@density3;

while($numelements3 < $numelements1){
	push(@density3,0);
	$numelements3=@density3;
}

$i = 0;
while(my $line = <INPUTFILE4>){
    chomp $line;
    my @linearray = split(" ", $line);
	push(@z4, $linearray[0]);
	push(@density4, $linearray[1]);
	if ($z1[$i] != $z4[$i]){
		unshift(@z4,$z1[$i]); #add element to TOP of array
		unshift(@density4,0);
	}	
	$i++;
}
my $numelements4=@z4;

while($numelements4 < $numelements1){
	push(@density4,0);
	$numelements4=@density4;
}

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 \"All\" \n";
print OUTPUTFILE "\@ s1 legend \"Water\" \n";
print OUTPUTFILE "\@ s2 legend \"POPC\" \n";
print OUTPUTFILE "\@ s3 legend \"Protein\" \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] \t $density3[$j] \t $density4[$j] \n";
}

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


