#!/usr/bin/perl

####### Edit here ######
  $prop_id = "o07116";               # Proposal ID
  $STARNAME = "OBJECT01_NL_1";       # Name of a set of integrations
# $STARNAME = "REF01_NL_1";          # Name of a set of integrations
  @N=(80286,80288,80290,80292,80294);# List of number 
# @N=(80370,80372);                  # List of number 
  $ymin = 98 - 5 ; $ymax = 98 + 5 ;  # y-position range 

  $RATIO_STARNAME = "REF01_NL_1";    # rationing star's STARNAME
# $RATIO_STARNAME = "none" ;         # none for no-rationing

  $gnuplot_command= "set yrange[0:1]"; # Gnuplot command what you like.
########################

# check and make directory
 if(! -e $STARNAME){
   $c=sprintf "mkdir $STARNAME";printf "$c\n"; system($c);
 }

# Find lambda relation
  $c = sprintf "q_sky_nlow /data/$prop_id/COMA%08d.fits 1 30-300 %d 1 2 |",$N[0],($ymax+$ymin)/2 ;
  open g,$c; $_=<g>; close g;
  if(/Best/){ printf "Cannot find x-lambda relation\n"; $A=1;$B=0; }
  else { chop; split ; $A=$_[5] ; $B = $_[7] ; }
  printf "lambda = $A * X + $B\n",$A,$B;

# Add all files
 if($#N==0){
   $c = sprintf "cp /data/$prop_id/COMQ%08d.fits $STARNAME/wrk1.fits",$N[0];
   printf "$c\n";system($c);
 }else{
   $c = sprintf "q_arith /data/$prop_id/COMQ%08d.fits + /data/$prop_id/COMQ%08d.fits $STARNAME/wrk1.fits",$N[0],$N[1];
   printf "$c\n";system($c);
 }
 for($i=2;$i<=$#N;$i++){
   $c=sprintf "q_arith $STARNAME/wrk1.fits + /data/$prop_id/COMQ%08d.fits $STARNAME/wrk2.fits",$N[$i];
   printf "$c\n";system($c);
   $c="cp $STARNAME/wrk2.fits $STARNAME/wrk1.fits";
   printf "$c\n";system($c);
 }

# List spectroscopic data
 $i=0;
 open h,"> $STARNAME/list.txt";
 $c = "q_list_stat $STARNAME/wrk1.fits 1 - $ymin:$ymax 1 | ";
 open g, $c;
 while(<g>){
   chop ; split ; $x = $_[1] ; $v = $_[4] ;
   $r = $A*$x + $B ;
   printf h "$r $v\n";
   $R[$i] = $r ; $V[$i] = $v ; $i++;
 }
 close g;
 close h;

# read rationing star
 if( -e "$RATIO_STARNAME/list.txt" ){
   $i=0;
   open h,"< $RATIO_STARNAME/list.txt" ;
   while(<h>){
     chop;split;$VR[$i]=$_[1];$i++;
   }
   close h;
 }else{
   for($i=0;$i<=$#V;$i++){$VR[$i]=1;}
 }

# plot data
 open g,"> tmp";
 for($i=0;$i<=$#V;$i++){
   if( $VR[$i] != 0 ){
     printf g "$R[$i] %e\n", $V[$i]/$VR[$i];
   }else {
     printf g "$R[$i] 0\n";
   }
 }
 close g;
 if($A!=1){
   $xrange = "xrange[7.5:13.5]";
 }else {
   $xrange = "xrange[1:320]";
 }

  open GNUPLOT, "| gnuplot -persist -geometry 640x480";
  print GNUPLOT <<gnuplot_Commands_DoneA;
  set title \"$STARNAME\"
  set xlabel \"Wavelength\"
  set $xrange
  $gnuplot_command
  set nokey
  plot \"tmp\" using 1:2 with lines 
gnuplot_Commands_DoneA