package cal_pipe2;

require Exporter;

@ISA      = qw(Exporter);
@EXPORT   = qw(specific_energy);

use strict;

sub specific_energy {

  my ($file, $outfile, %head_vals) = @_;

  my $ir_state = $head_vals{ir_power_state} eq "ON";
  my $vis_state = $head_vals{vis_power_state} eq "ON";

  my ($bands, $this_function, $cmd);

  $bands = "1-352";
  my $temptbl = "temp_table";
  my $temptbl2 = "temp_table.tbl";

  my $cube_mode = $head_vals{instrument_mode_id};

  my  $ir_perf  = "$ENV{caldir}/calib/ir_perf.qub";
  my  $vis_perf = "$ENV{caldir}/calib/vis_perf.qub";
  my  $wavecal  = "$ENV{caldir}/calib/wavecal.qub";

  if ($head_vals{ssf} eq "ON") {
    $ir_perf  = "$ENV{caldir}/calib/ir_perf_ss.qub";
    $vis_perf = "$ENV{caldir}/calib/vis_perf_ss.qub";
    $wavecal  = "$ENV{caldir}/calib/wavecal_ss.qub";
  }


  # Paste an extra line on the bottom of the cube so we can use multcs
  # without getting the multiplier spectrum stuck in the cube.
  $cmd = "cubestack from1=$file" .
     "  sfrom1='\"::\"'" .
     "  from2=$file" .
     "  sfrom2='\":1:\"'" .
     "  to=$outfile" .
     "  otype=3" .
     "  orange=(0.0,0.0)" .
     "  HOPT=1";
  print ("$cmd\n");
  system($cmd);


  # Multiply IR cube by IR performance model.
  if ($ir_state) {
    $cmd = "cubespec from=$ir_perf" .
	   "  sfrom='\"1:1:$bands\"'" .
	   "  tblto=$temptbl";
    print ("$cmd\n");
    system($cmd);

    system("mv $outfile temp.cub");
    $cmd = "multcs from=temp.cub" .
	   "  sfrom='\"::\"'" .
	   "  to=$outfile" .
	   "  otype=3" .
	   "  orange='\"-10000.0,100000000.0\"'" .
	   "  tblfrom=$temptbl";
    print ("$cmd\n");
    system($cmd);
    system("rm $temptbl2");
    system("rm temp.cub");

    if ($vis_state) {
      $this_function = "'(\"F1*(BAND<97)+F1/2*(BAND>=97)\")'";
    } else {
      $this_function = "'(\"F1/2\")'";
    }

    # If the gain is HIGH, divide by 2.
    if ($head_vals{ir_gain_mode} eq "HIGH") {
      system("mv $outfile temp.cub");
      $cmd = "fx from=(temp.cub)" .
	     " sfrom='(\"::\")'" .
	     " function=$this_function" .
	     " to=$outfile otype=3";
      system($cmd);
      system("rm temp.cub");
    }
  }

  system("rm isis_param.out");

  # Multiply VIS cube by VIS performance model.
  if ($vis_state) {
    system("mv $outfile temp.cub");
    $cmd = "cubespec from=$vis_perf" .
	   "  sfrom='\"1:1:$bands\"'" .
	   "  tblto=$temptbl";
    print ("$cmd\n");
    system($cmd);

    $cmd = "multcs from=(temp.cub)" .
	   "  sfrom='\"::\"'" .
	   "  to=$outfile" .
	   "  otype=3" .
	   "  orange='\"-10000.0,100000000.0\"'" .
	   "  tblfrom=$temptbl";
    print ("$cmd\n");
    system($cmd);
    system("rm $temptbl2");
    system("rm temp.cub");

    if ($ir_state) {
      $this_function = "'(\"F1/2*(BAND<97)+F1*(BAND>=97)\")'";
    } else {
      $this_function = "'(\"F1/2\")'";
    }

    # If the gain is HIGH, divide by 2.
    if ($head_vals{vis_gain_mode} eq "HIGH") {
      system("mv $outfile temp.cub");
      $cmd = "fx from=(temp.cub)" .
	     " sfrom='(\"::\")'" .
	     " function=$this_function" .
	     " to=$outfile otype=3";
      system($cmd);
      system("rm temp.cub");
    }

    my $v_high_res_mult = 3.0;

    if ($ir_state) {
      $this_function = "'(\"F1*(BAND<97)*$v_high_res_mult+F1*(BAND>=97)\")'";
    } else {
      $this_function = "'(\"F1*$v_high_res_mult\")'";
    }

    if ($head_vals{vis_sampling_mode} eq "HI-RES") {
      system("mv $outfile temp.cub");
      $cmd = "fx from=(temp.cub)" .
	     " sfrom='(\"::\")'" .
	     " function=$this_function" .
	     " propfile=1" .
	     " to=$outfile otype=3";
      system($cmd);
      system("rm temp.cub");
    }
  }

  system("rm isis_param.out");

  # Multiply the data by (h*c)/(lambda*A*omega*delta_lamda) to get
  # specific enengy.
  system("mv $outfile temp.cub");
  $cmd = "cubespec from=$wavecal" .
	 "  sfrom='\"1:1:$bands\"'" .
	 "  tblto=$temptbl";
  print ("$cmd\n");
  system($cmd);

  $cmd = "multcs from=(temp.cub)" .
	 "  sfrom='\"::\"'" .
	 "  to=$outfile" .
	 "  otype=3" .
	 "  orange='\"-10000.0,100000000.0\"'" .
	 "  tblfrom=$temptbl";
  print ("$cmd\n");
  system($cmd);
  system("rm $temptbl2");
  system("rm temp.cub");

  system("rm isis_param.out");

  # Divide the vis part of the cube by the vis exposure (in seconds)
  if ($vis_state) {

    my $vis_exp = $head_vals{vis_exp}/1000.0;

    if ($ir_state) {
      $this_function = "'(\"F1*(BAND<97)/$vis_exp+F1*(BAND>=97)\")'";
    } else {
      $this_function = "'(\"F1/$vis_exp\")'";
    }

    system("mv $outfile temp.cub");
    $cmd = "fx from=(temp.cub)" .
	   " sfrom='(\"::\")'" .
	   " function=$this_function" .
	   " propfile=1" .
	   " to=$outfile otype=3";
    system($cmd);
    system("rm temp.cub");
  }

  # Divide the ir part of the cube by the ir exposure (in seconds)
  if ($ir_state) {
    my $ir_exp = $head_vals{ir_exp}/1000.0;
    $ir_exp = $ir_exp - .004;

    if ($vis_state) {
      $this_function = "'(\"F1*(BAND<97)+F1/$ir_exp*(BAND>=97)\")'";
    } else {
      $this_function = "'(\"F1/$ir_exp\")'";
    }

    system("mv $outfile temp.cub");
    $cmd = "fx from=(temp.cub)" .
	   " sfrom='(\"::\")'" .
	   " function=$this_function" .
	   " propfile=1" .
	   " to=$outfile otype=3";
    system($cmd);
    system("rm temp.cub");
  }

  # Now we have ergs/second.

  $cmd = "labels from=$outfile" .
         "  list=NONE" .
	 "  dest=FILE" .
         "  listfile=listfile" .
	 "  grpnam=CALIBRATION" .
         "  keyword=SPECIFIC_ENERGY" .
         "  keyvalue=V1.0" .
         "  keytype=STRING";
  system($cmd);
  system("rm listfile.lis");

  # Trim off that extra line we added earlier.
  system("mv $outfile temp.cub");
  my $lns = $head_vals{'lines'};
  $cmd = "transpose from=temp.cub sfrom='(\":1-$lns:\")' to=$outfile " .
          "order=BSQ\n";
  print("$cmd");
  system($cmd);
  system("rm temp.cub");

  system("rm isis_param.out");
}