#!/usr/bin/env perl

use POSIX;

$default_linewidth=3;  # in km/s
$cutoff=0.25;             # limiting sensitivity in K km/s
$initdens=3.5;         # starting density of log10(H2).
$Tscale=1.5;
$peakCO = '12co_peak.csv';   # name of 12CO mom=-2 map
$meanCO = '12co_mean.csv';   # name of 12CO mom=0 map
$mean13CO = '13co_mean.csv';   # name of 13CO mom=0 map
$widthCO = '12co_width.csv'; # name of 12CO mom=2 map
$width13CO = '13co_width.csv'; # name of 13CO mom=2 map

$moleID1 = '../grids/co2-1';      # name of the molecular grids
$moleID2 = '../grids/13co2-1';    

# You shouldn't need to change anything below this line

@colnames = ("1e15", "2e15", "5e15", "1e16", "2e16", "5e16", "1e17", "2e17", "5e17", "1e18");
local $| = 1;
$dens=$initdens;
$i=0;

# first read peak line brightness to get Tkin
print "Reading 12CO for Tkin grid... \n";
open(FILE,$peakCO) ||  die "Data file not found";
@dummy=<FILE>;
close(FILE);
&Remap;
@Tkin = @array;
$i=0;
foreach (@Tkin)
{
    $Tkin[$i]*=$Tscale;
    if($Tkin[$i] < 10)
    {
        $Tkin[$i]=10;
    }
    $i++;
}

print "Reading integrated intensities and widths for 12CO...\n";

open(FILE,$meanCO) ||  die "Data file not found";
@dummy=<FILE>;
close(FILE);
&Remap;
@intensity12 = @array;

open(FILE,$widthCO) ||  die "Data file not found";
@dummy=<FILE>;
close(FILE);
&Remap;
@width12=@array;
$i=0;

print "Reading integrated intensities and widths for 13CO...\n";

open(FILE,$mean13CO) ||  die "Data file not found";
@dummy=<FILE>;
close(FILE);
&Remap;
@intensity13 = @array;

open(FILE,$width13CO) ||  die "Data file not found";
@dummy=<FILE>;
close(FILE);
&Remap;
@width13=@array;
$i=0;

print "Reading Escape Probability grids...\n";
    $j=0;
    foreach $col (@colnames) 
    {
        $filename = $moleID1."_".$col.".dat";
        open(GRID, $filename) ||  die "grid file not found";
        @{"grid1".$j}=<GRID>;
        close(GRID);
        $j++;
    }
    $jmax=$j;

    $j=0;
    foreach $col (@colnames) 
    {
        $filename = $moleID2."_".$col.".dat";
        open(GRID, $filename) ||  die "grid file not found";
        @{"grid2".$j}=<GRID>;
        close(GRID);
        $j++;
    }


print "Now solving for column density...\nProgress (", $max, " total):";

open(FILE, "> output.csv");
open(FILE2, "> radiation.csv");

for($i=0;$i<=$max;$i++)
{
  if($intensity13[$i] < $cutoff)
  {
      if($intensity12[$i] < $cutoff)
      {
          $N = 0;
          print FILE $N, "\n";
          print FILE2 "0\n";
      }
      else
      {
          $temp=$Tkin[$i]*2.0;
          $n=$dens-0.5;
#	  $n=$dens;
          $intensity=$intensity12[$i];
          $grid="grid1";
          &Solver;
          print FILE $N, "\n";
          print FILE2 pow(($Tkin[$i]/10.0),5), "\n";
      }
  }
  else
  {
      $temp=$Tkin[$i];
      $n=$dens;
      $intensity=$intensity13[$i];
      $grid="grid2";
      &Solver;
      $N*=50;
      print FILE $N, "\n";  
      print FILE2 pow(($Tkin[$i]/10.0),5), "\n";
  }
  if($i % 1000 == 0)
  {
      print " ", int($i/1000);
  }
}
close(FILE2);
close(FILE);
print "\n";


sub Remap
{
    $max=0;
    $i=0;
    foreach (@dummy)
    {
        @line = split();
        $j=0;
        foreach (@line)
        {
            $array[8*$i+$j]=$_;
            $max++;
            $j++;
        }
        $i++;
    }
}


sub Solver
{
    $found=0;

    # optional:  scale density with 13CO column
#    $scaling = $intensity13[$i];
#    if($scaling <= 1.0)
#    {
#        $scaling=1.0;
#    }
#    $scaling=log10($scaling);
#    $n += $scaling;

    # find best linewidth at this position
    if($intensity13[$i] < $cutoff*2.0)
    {
        if($intensity12[$i] < $cutoff*2.0) 
        {
            $width=$default_linewidth;
            $jmax=10;
        }
        else
        {
            $width=@width12[$i];
            $jmax=10;
        }        
    }
    else 
    {
        $width=@width13[$i];
        $jmax=8;
    }
    if($width < 0.1)
    {
        $width=$default_linewidth;
    }
  OUTER: for($j=0; $j < $jmax; $j++)
  {
      $count=0;
      foreach $buf (@{$grid.$j})
      {
          @line=split(' ', $buf);
          if(($line[0] >= $temp) && ($found < 1))
          {
              if(($line[1] >= $n) && ($count < 1))
              {
                  $count++;
                  #print "\t", $line[2], " ", $intensity, " ",$j, "\n";
                  if($line[2] > $intensity/$width)
                  {
                      $found++;
                      $ratio = $intensity/$width/$line[2];
                      $N = $ratio*$colnames[$j]*$width;
                      last OUTER;
                  }
              }
          }
      }
  }
    if($found < 1)
    {
        $N=0;
    }
#    $n=$initdens;
}
