#!/usr/bin/perl

  $prop_id = "o07116";
  $STARNAME = "TEST_A";
  @N=(80247);

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

# loop of file 
 for($ii=0;$ii<=$#N;$ii++){
# STEP0 copy file
  $coma = sprintf "COMA000%d.fits",$N[$ii];
  $comq = sprintf "COMQ000%d.fits",$N[$ii];
  if( ! -e $STARNAME."/".$coma ){
   $c=sprintf "cp /data/$prop_id/$coma $STARNAME";system($c); 
  }
  if( ! -e $STARNAME."/".$comq ){ 
   $c=sprintf "cp /data/$prop_id/$comq $STARNAME";system($c); 
  }
 
# STEP1 find min-max position
  $xmin = 20 ; $xmax = 300 ; $ymin = 20 ; $ymax = 220 ;
  $c=sprintf "q_minmax $STARNAME/$comq 1 $xmin:$xmax $ymin:$ymax 1 -p |"; 
  open g,$c ; $_=<g> ; close g; 
  if(/NAXIS2/){
   $c=sprintf "q_minmax $STARNAME/$comq 1 : : 1 -p |";
   open g,$c ; $_=<g> ; close g;
  }
  split; $minx = $_[3]; $miny = $_[4]; $maxx = $_[10]; $maxy = $_[11];

  printf "Star at $maxx,$maxy and $minx,$miny\n" ;

# STEP2 find center position
  Q_PHOTO($STARNAME."/".$comq,$maxx-15,$maxx+15,$maxy-15,$maxy+15,
              \$maxxc,\$maxyc,\$sn,\$fwhm) ;
  printf "PLUS-side  Star at $maxxc,$maxyc  SN= $sn\n";
  if($sn<3 && $sn >-3 ){ print "DIE\n";die "DIE : Too low SN $sn\n"; }
  Q_PHOTO($STARNAME."/".$comq,$minx-15,$minx+15,$miny-15,$miny+15,
              \$minxc,\$minyc,\$sn,\$fwhm) ;
  printf "MINUS-side Star at $minxc,$minyc  SN= $sn\n";
  if($sn<3 && $sn >-3 ){ print "DIE\n";die "DIE : Too low SN $sn\n"; }
  
# STEP3 find out the region
  for($side=0;$side<2;++$side){
    if($side==0){$ss="PLUS";$cx=$maxxc;$cy=$maxyc;}
    else{$ss="MINUS";$cx=$minxc;$cy=$minyc;}
    $x0[$side]=sprintf"%d",$cx; $y0[$side]=sprintf"%d",$cy;
    $x0[$side]=$x0[$side]-32;$x1[$side]=$x0[$side]+64;
    $y0[$side]=$y0[$side]-32;$y1[$side]=$y0[$side]+64; 
    if($x1[$side]>320){$x1[$side]=320;$x0[$side]=320-64;}
    if($x0[$side]<1){$x0[$side]=1;$x1[$side]=65;}
    if($y1[$side]>240){$y1[$side]=240;$y0[$side]=240-64;}
    if($y0[$side]<1){$y0[$side]=1;$y1[$side]=65;}
    printf "$ss cut $x0[$side]-$x1[$side] $y0[$side]-$y1[$side]\n";
  }

# STEP4 Check the mode
  $mode = Q_HEAD($STARNAME."/".$coma,"Q_CHAM");
  $nexp = Q_HEAD($STARNAME."/".$coma,"Q_CHEB"); # exposure / beam
  $nchop= Q_HEAD($STARNAME."/".$coma,"Q_CHCN"); # chop / frame
  $nz   = Q_HEAD($STARNAME."/".$coma,"NAXIS3");
  printf "MODE= $mode NEXP= $nexp NCHOP= $nchop NZ= $nz\n";
  if($mode==0 && $nz == $nexp*$nchop){
    $nnexp=$nexp; printf "ROW mode OK\n" ;
  }elsif($mode==1 && $nz == $nchop ){
    $nnexp=1; printf "ADD mode OK\n" ;
  }else{ print "DIE\n";die "DIE : Unsupported mode\n";}

# STEP5 Cut sub image
  for($side=0;$side<2;++$side){ 
    if($side==0){$f="$STARNAME/$coma.P";}
    else        {$f="$STARNAME/$coma.M";}
    if(! -e $f){
      $c = sprintf "q_list_stat $STARNAME/$coma 1 %d-%d %d-%d - $f",
          $x0[$side],$x1[$side],$y0[$side],$y1[$side];
      system($c);
    }
  }

# STEP6 shift and add
  if(! -e "$STARNAME/$coma.sa") {
    $kk=0;
    open h, "> $STARNAME/$coma.sa";
    open g1, "> wrk1";
    open g2, "> wrk2";
    for($j=0;$j<$nchop;$j=$j+2){
      $c=sprintf "q_list_stat $STARNAME/$coma.P 1 - - %d:%d tmppb",
               ($j+1)*$nnexp+1,($j+1)*$nnexp+$nnexp;
      system($c);
      $c=sprintf "q_list_stat $STARNAME/$coma.M 1 - - %d:%d tmpmb",
               ($j)*$nnexp+1,($j)*$nnexp+$nnexp;
      system($c);
      for($k=0;$k<$nnexp;++$k){
        $c=sprintf"q_list_stat $STARNAME/$coma.P 1 - - %d tmppa",
             $j*$nnexp+$k+1;
        system($c);
        $c=sprintf"q_list_stat $STARNAME/$coma.M 1 - - %d tmpma",
             ($j+1)*$nnexp+$k+1;
        system($c);
        $c=sprintf "q_arith tmppa - tmppb tmppc";system($c);
        $c=sprintf "q_arith tmpma - tmpmb tmpmc";system($c);
        Q_PHOTO("tmppc",18,48,18,48,\$pcx,\$pcy,\$psn,\$pfwhm) ;
        if($psn<1.1){print "DIE\n";die "DIE : Too low SN<1.1 $psn\n"; } 
        printf  "%03d PLUS  chop= %02d exp= %02d SN= %6.2f FWHM= %3.1f center $pcx $pcy\n",$j*$nnexp+$k,$j,$k,$psn,$pfwhm;
        printf h "%03d PLUS  chop= %02d exp= %02d SN= %6.2f FWHM= %3.1f center $pcx $pcy\n",$j*$nnexp+$k,$j,$k,$psn,$pfwhm;
        Q_PHOTO("tmpmc",18,48,18,48,\$mcx,\$mcy,\$msn,\$mfwhm) ;
        if($msn<1.1){print "DIE\n";die "DIE : Too low SN<1.1 $msn\n"; }
        printf  "%03d MINUS chop= %02d exp= %02d SN= %6.2f FWHM= %3.1f center $mcx $mcy\n",($j+1)*$nnexp+$k,$j,$k,$msn,$mfwhm;
        printf h "%03d MINUS chop= %02d exp= %02d SN= %6.2f FWHM= %3.1f center $mcx $mcy\n",($j+1)*$nnexp+$k,$j,$k,$msn,$mfwhm;
        $c=sprintf"q_expand tmppc 10 tmppd";system($c);
        $c=sprintf"q_expand tmpmc 10 tmpmd";system($c);
        $c=sprintf"q_list_stat tmppd 1 %d-%d %d-%d 1 tmppe.%03d",
              $pcx*10-120,$pcx*10+120,$pcy*10-120,$pcy*10+120,$kk ;
        system($c);
        $c=sprintf"q_list_stat tmpmd 1 %d-%d %d-%d 1 tmpme.%03d",
              $mcx*10-120,$mcx*10+120,$mcy*10-120,$mcy*10+120,$kk ;
        system($c);
        printf g1 "tmppe.%03d\n",$kk; printf g2 "tmpme.%03d\n",$kk;
        $kk++;
      }
    }
    close g1 ; close g2 ; close h ;
    $c=sprintf"q_fcombine \@wrk1 ave=$STARNAME/$coma.P.sa"; 
    printf "$c\n";system($c);
    $c=sprintf"q_fcombine \@wrk2 ave=$STARNAME/$coma.M.sa"; 
    printf "$c\n";system($c);
    system("/usr/bin/rm tmp* wrk*");
  }

}
exit ;

#################################
sub Q_PHOTO{
  my ($file, $x0, $x1, $y0, $y1, $xc,$yc,$sn,$fwhm) = @_ ;
  my $c=sprintf "q_photo_gif $file 1 $x0-$x1 $y0-$y1 1 2>/dev/null |" ;
  my ($a);
# printf "$c\n";
  open g,$c ;
  while(<g>){#print ;
    if(/FWHM/){ 
      s/=/ /g;s/\(/ /;s/\)/ /;split; $$xc = $_[3],$$yc = $_[4] ; 
      $a=$_[1];$$fwhm=$_[6];
    }
    if(/^SN/){chop;s/=/ /;split;$$sn=$_[1];}
  } 
  close g;
  if(-e "qplot.gif") { system("/usr/bin/rm qplot.gif"); }
}
#################################
sub Q_HEAD {
  my $key=pop; my $file=pop;
  open q_head,"q_head $file $key |";$_=<q_head>;close q_head; chop;split;
  $_[0];
}