#!/usr/bin/perl -w

# Simple script to convert wave analysis "strain" output to usable
# hplus, hcross for data analysis purposes
# Adapted from c++ code of Sean McWilliams

# This code asusmes that only the (2,2) and (2,-2) modes contribute
# to the detected strain, and moreover, that they contribute equally.

use FileHandle;
use Math::Trig;

if(@ARGV != 5)
{
  print "Usage: $0 <input data> <total redshifted mass (M_Sun)> <luminosity distance (Mpc)> <inclination (degs)>
<polarization (degs)> \n";
  exit;
}

$MADM = 1.0;
$PI = 2.0*acos(0.0);
#$C = 2.9979e8;	#speed of light         [m/s]
#$G = 6.6732e-11; #Newton's gravitational constant	[m^3/kg/s^2]
$t_conv = 4.9169e-6; #convert t = [M] to t = [sec] 
$h_conv = 4.7804e-20; #convert h = [M]/[r] to h dimensionless

$infile = $ARGV[0];
$M = $ARGV[1];
$D = $ARGV[2];
$iota = $PI*$ARGV[3]/180.0;
$pol = $PI*$ARGV[4]/180.0;

print "M = ",$M,", D = ",$D,", iota = ",$iota,", pol = ",$pol,"\n";

@raw_time = ();
@hRe = ();
@hIm = ();

$fileroot = substr($infile, 0, -4);
#$outparams = $fileroot.".param";
$outfile = $fileroot."_at_detector.dat";


open(IN_FILE,   "<$infile") || die "Unable to find file \"$infile\".";

$mode2pol = sqrt(5.0/(4.0*$PI))*$h_conv*$M/$D/$MADM;

$ii = 0;
while (<IN_FILE>) {
  $line = $_;
    
  if($line eq "\n") { # If line is empty line
    # dump
  } elsif(/^\#/) { # If line has comments
    # dump also
  } else {

    chop($line);
    @data = split(/\s/,$line);
    $raw_time[$ii] = $data[0];
    $hRe[$ii] = $data[1];
    $hIm[$ii] = $data[2];

    $ii++; 
  }
    
}

close(IN_FILE);

$dt = $raw_time[1] - $raw_time[0];
$dt = $t_conv*$M*$dt/$MADM;

#@time_val = ();
@hplus_val = ();
@hcross_val = ();

open(OUT_FILE,   ">$outfile") || die "Unable to create file \"$outfile\".";

print(OUT_FILE "#time\th plus\th cross\n");

for $jj (0 .. ($ii-1)) {

#  $time_val[$jj] = $t_conv*$M*$raw_time[$jj]/$MADM;

  $tnow = $jj*$dt;
  $hplus_val[$jj] = $mode2pol*(cos(2*$pol)*$hRe[$jj]+sin(2*$pol)*$hIm[$jj])*(1 + (cos($iota))*(cos($iota)))/2;
  $hcross_val[$jj] = $mode2pol*(-sin(2*$pol)*$hRe[$jj]+cos(2*$pol)*$hIm[$jj])*cos($iota);

# print(OUT_FILE $tnow,"        ",$hplus_val[$jj],"        ",$hcross_val[$jj],"\n");
  printf(OUT_FILE "%.8e\t%.8e\t%.8e\n",$tnow,$hplus_val[$jj],$hcross_val[$jj]);
}

close(OUT_FILE);