#!/usr/bin/perl 
use Math::Trig;

$PI = 3.1415926535;
$MAGDIFF = 1;

#################### Example of galfit.pl control parameters  #################

# $ACSIMG = "s1z09.fits"  # The ACS image name to fit.
# $SEXCAT = "s1z09.cat"   # The Sextractor image catalogue name.
# $PSFLOC = "psfs"        # The directory storing all the PSFs.
# $STARTNUM = 1 	  # Object starting number to fit.
# $STOPNUM = 999 	  # Object stopping number to fit.
# $CONSFILE = "constraints"  # The constraint file name.
# $TEMPLOC = "galfit-tempfits"  # The directory storing image sub-tiles.
# $FITFUNC = "sersic"     # Fit a single sersic or do B/D decomposition?
# $REGION = "A"           # (A) Whole image catalog or, (B) region of the image;
# $BOUNDARY = "  "        # Boundary of the object region if (B);
# $SKYVAL = 0.	 	  # Sky value added back to an image that SExtractor knows nothing about.
# $MRANGE = "0. 99."      # Acceptable magnitude range;
# $STARGAL = "0.0 1.0"    # Range of good galaxies: 0.0 (galaxy) to 1.0 (star);
# $FWHMPSF = 0.05         # FWHM of the PSF for cosmic ray rejection.
# $MAGZPT = 23.64         # Photometric zeropoint.
# $PLATE = 0.03           # Plate scale
# $NCOL, $NROW = "512 512"  # The image size of the sub-tiles.
# $TNAME = "temp"	  # The prefix name of the image sub-tiles.
# $BUFFER = "200"	  # The padding size around the image edges

############# Read in the control parameters as defined above  ###############

$INFILE = @ARGV[0];

open (IN1, $INFILE) || die "Can't open $INFILE $!";

$line=<IN1>;
($ACSIMG) = split ' ', $line;	  # The ACS image name to fit.
$line=<IN1>;
($SEXCAT) = split ' ', $line;	  # The sextractor catalog.
$line=<IN1>;
($PSFLOC) = split ' ', $line;	  # The directory storing all the PSFs.
$line=<IN1>;
($STARTNUM) = split ' ', $line;   # Object starting number to fit.
$line=<IN1>;
($STOPNUM) = split ' ', $line;    # Object stopping number to fit.
$line=<IN1>;
($CONSFILE) = split ' ', $line;	  # The constraint file name.
$line=<IN1>;
($TEMPLOC) = split ' ', $line;    # The directory storing image sub-tiles.
$line=<IN1>;
($FITFUNC) = split ' ', $line;    # Fit a single sersic or do B/D decomposition?
$line=<IN1>;
($REGION) = split ' ', $line;     # (A) Whole image catalog or, (B) region of an image.
$line=<IN1>;
($BOUNDARY) = split ' ', $line;   # Boundary of the object region if (B).
$line=<IN1>;
($SKYVAL) = split ' ', $line;     # Sky value added back to an image that SExtractor knows nothing about.
$line=<IN1>;
($MRANGE) = split '#', $line;     # Acceptable magnitude range.
$line=<IN1>;
($STARGAL) = split '#', $line;    # Range of good galaxies: 0.0 (galaxy) to 1.0 (star).
$line=<IN1>;
($FWHMPSF) = split ' ', $line;    # FWHM of the PSF for cosmic ray rejection.
$line=<IN1>;
($MAGZPT) = split ' ', $line;     # Photometric zeropoint.
$line=<IN1>;
($PLATE) = split ' ', $line;      # The image plate scale.
$line=<IN1>;
($NCOL, $NROW) = split ' ', $line;  # The image size of the sub-tiles.
$line=<IN1>;
($TNAME) = split ' ', $line;      # The prefix name of the image sub-tiles.
$line=<IN1>;
($NSPLIT) = split ' ', $line;     # The number of subtiles along one axis.
$line=<IN1>;
($BUFFER) = split ' ', $line;	  # The padding size around the image edges.

close (IN1);

if ($REGION eq "B") { ($x1,$y1,$x2,$y2)=split " ",$BOUNDARY; }
($brtcut,$ftcut)=split " ",$MRANGE;
($cut1,$cut2)=split " ",$STARGAL;
$minfwhm = (0.95*$FWHMPSF)/0.05;      #rids Crays
$xpansz = int ($NCOL / $NSPLIT);
$ypansz = int ($NROW / $NSPLIT);

##############################################################################

print "=============================================================\n";
print " This script takes the SExtractor output, generated by Dan   \n";
print " McIntosh and formats it for GALFIT.  The sextractor output  \n";
print " has 17 columns.						    \n";
print "=============================================================\n";

open (IN2, $SEXCAT) || die "Can't open $SEXCAT $!";

###############  Read in sextractor data, weed out bad objs  #################

$i=1;            #good array element
$ii=0;           #overall counting


open (OUT2, "> $TEMPLOC/offsets");

while($line=<IN2>){
  chomp($line);
  ($num,$f,$fe,$mg,$mgerr,$rk,$bkgd,$ia,$xx,$yy,$ra,$dec,$pa,$ell,$fw,
						$flg,$idx)=split ' ',$line;
  if ($num ne "\#") {
    $ii++;
    # std. flag and cosmic ray culling
    if ($flg < 4 && $fw > $minfwhm) { 
#     if ($mg <= $ftcut && $mg >= $brtcut) {
        if ($idx <= $cut2 && $idx >= $cut1) {
          if ($REGION eq "A" || ($REGION eq "B" && $xx >= $x1 && $xx <= $x2 && $yy >= $y1 && $yy <= $y2)) {

	    $objid[$i] = $num;
            $flux[$i]=$f; $fluxerr[$i]=$fe;
            $mag[$i]=$mg; $magerr[$i]=$mgerr;
            $Re[$i]=0.3*sqrt($ia);        # fudge
            $bkg[$i]=$bkgd;$isoa[$i]=$ia;
            $x[$i]=$xx; $y[$i]=$yy;
            $alpha[$i]=$ra; $delta[$i]=$dec;
            $fwhm[$i]=$fw;
            #elong[$i]=$elo; 
            $axrat[$i]=1.-$ell;
            $flag[$i]=$flg; $index[$i]=$idx;

       # Sextractor to galfit orientation
	    $theta[$i] = $pa - 90;
            $i++;
          } #img. boundary
        } #star/gal
#     } #mags
    } #flag,crays
  }
}
$totobj=$i - 1;

open (OUT3, "> galfit.crashes");
open (OUT4, "> galfit.fitted");
run_galfit ();

print "Number of good sources in $SEXCAT is $totobj from total of $ii\n";

close (OUT4);
close (OUT3);
close (OUT2);
close (IN2) || die "Can't close $SEXCAT $!";


############  New subroutine to create GALFIT input files  ###################


sub run_galfit {

    my ($i, $j, $ix, $iy, $ncomp, $xsize, $ysize, $radck);

    my $NRe = 2.;

    $i = 0;
    while ($i <= $totobj && $errno != 256) {
	$i++;
        if ($objid[$i] >= $STARTNUM && $objid[$i] <= $STOPNUM &&
                            ($mag[$i] <= $ftcut && $mag[$i] >= $brtcut)) {
	    $sobj = $objid[$i];
	    $ncomp = 1;
            $parfile = "obj" . $objid[$i];
            open (OUT1, "> $parfile");

         # Experimental image size
	    $majx = abs(sqrt($isoa[$i]) * 4. * sin ($theta[$i]/180.*$PI));
	    $majy = abs(sqrt($isoa[$i]) * 4. * cos ($theta[$i]/180.*$PI));
	    $minx = $axrat[$i]* sqrt($isoa[$i]) * 4. * 
				abs(cos (($theta[$i]+$PI/2)/180.*$PI));
	    $miny = $axrat[$i]* sqrt($isoa[$i]) * 4. * 
				abs(sin (($theta[$i]+$PI/2)/180.*$PI));

	    $xsize = $majx;
	    if ($minx >= $majx) {
	        $xsize = $minx;
	    };
							
	    if ($xsize <= 30) {
		$xsize = 30;
            };

	    $ysize = $majy;
	    if ($miny >= $majy) {
	        $ysize = $miny;
	    };
							
	    if ($ysize <= 30) {
		$ysize = 30;
	    };

         # Figure out which subpanel this current object falls into
	    $ix = int ($x[$i] / $xpansz) + 1;
	    $iy = int ($y[$i] / $ypansz) + 1;
	    $fitname = $TNAME . "-" . $ix . "-" . $iy . ".fits";
	    $rmsname = "rms-" . $ix . "-" . $iy . ".fits";

         # Now make sure the object coordinate in the subpanel is correct
            if ($ix == 1) {
                $xbuffer = 0;
            } else {
                $xbuffer = $BUFFER;
            };

            if ($iy == 1) {
                $ybuffer = 0;
            } else {
	        $ybuffer = $BUFFER;
            };

         # Figure out which PSF is closest to object center
            $psffile = psfdist ($x[$i], $y[$i]);

         # Construct postage stamp name
            ($ra_cel,$dec_cel)=deg2celestial($alpha[$i],$delta[$i]);
            $outname = uniqname ($ra_cel, $dec_cel);

            ####################################
            #  Create GALFIT input file header #
            ####################################

            print OUT1 ("================================================================================\n");
            print OUT1 ("# IMAGE PARAMETERS \n");
            print OUT1 ("A) $TEMPLOC/$fitname   # Input Data image (FITS file) \n");
            print OUT1 ("B) $outname-out.fits   # Output data image block \n");
            print OUT1 ("C) $TEMPLOC/$rmsname   # Noise image name (made from data if blank or \"none\") \n");
            print OUT1 ("D) $psffile            # Input PSF image for convolution (FITS file) \n");
            print OUT1 ("E) 1                   # PSF oversampling factor relative to data \n");
            print OUT1 ("F) none                # Bad pixel mask (FITS image or ASCII coord list) \n");
   	    print OUT1 ("G) $CONSFILE           # Parameter constraint file \n");
            print OUT1 ("J) $MAGZPT             # Magnitude photometric zeropoint \n");
            print OUT1 ("K) $PLATE $PLATE       # Plate scale (dx dy). \n");
            print OUT1 ("O) both                # Display type (regular, curses, both) \n");
            print OUT1 ("P) 0                   # Create ouput only? (1=yes; 0=optimize) \n");
            print OUT1 ("S) 0                   # Modify/create objects interactively? \n");
            print OUT1 (" \n");
            print OUT1 ("# INITIAL FITTING PARAMETERS \n");
            print OUT1 ("# \n");
            print OUT1 ("#   For object type, allowed functions are: sersic, nuker, \n");
            print OUT1 ("#                       expdisk, devauc, moffat, gaussian. \n");
            print OUT1 (" \n");
            print OUT1 ("# Objtype:      Fit?         Parameters \n");
            print OUT1 (" \n");

         # Calculate the (x,y) position of the current object relative to
         # the tile in which it lives.

            $offx = ($ix-1) * $xpansz - $xbuffer;
  	    $offy = ($iy-1) * $ypansz - $ybuffer;

	    $xfit = $x[$i] - $offx;
	    $yfit = $y[$i] - $offy;
	    $bkgnd = $bkg[$i] + $SKYVAL;

         # Calculate fitting box needed to plug into galfit header:

            $xmin = int ($xfit - $xsize);
            if ($xmin <= 0) {$xmin = 1};
            $xmax = int ($xfit + $xsize);
            $ymin = int ($yfit - $ysize);
            if ($ymin <= 0) {$ymin = 1};
            $ymax = int ($yfit + $ysize);

     # Figure out the nearest neighbors that needed to be deblended
     # together.

	    $ncomp = components ($ncomp, $i, $i, $offx, $offy);
            for ($j=1; $j <= $totobj; $j++) {
                if ($i != $j) {
		    $dx = $x[$i] - $x[$j];
		    $dy = $y[$i] - $y[$j];
		    $dist = sqrt($dx**2 + $dy**2);

     #  Fit all objects that fall within a certain radius of the central obj.:

		    $radck = sqrt($isoa[$j]);
		    $angle = atan2 ($dy, $dx);
		    $phi = $PI/2. - ($theta[$j]/180.*$PI + $PI/2. - $angle);
                    $thetap = atan2 ($axrat[$j] * tan ($phi), 1);

		    $isorad = $radck * sqrt (($axrat[$j] * cos ($thetap))**2 + 
						           sin($thetap)**2);

     # If object is within the fitted image region...
		    if ((($x[$j]-$offx >= $xmin) && ($x[$j]-$offx <= $xmax) && 
		         ($y[$j]-$offy >= $ymin) && ($y[$j]-$offy <= $ymax) &&
     # but is greater than 2 Re away, and is not more than 3 mag fainter, or
			 (($dist > 2 * $Re[$i] && $mag[$j]-3 <= $mag[$i]) ||
     # is less than 2 Re away and not more than 5 mag fainter, or
			  ($dist < 2 * $Re[$i] && $mag[$j]-5 <= $mag[$i]))) ||
     # is outside image, but is about 3 Re away and 1 mag brighter.
		        (($x[$j]-$offx >= $xmin-$isorad) && ($x[$j]-$offx <= $xmax+$isorad) &&
		         ($y[$j]-$offy >= $ymin-$isorad) && ($y[$j]-$offy <= $ymax+$isorad) &&
			                 $mag[$j] + $MAGDIFF <= $mag[$i]))  {

# print "$radck, $dist -- $theta[$j], $thetap, $phi -- $dy, $dx, $angle\n";


		        $ncomp++;
		        $ncomp = components ($ncomp, $i, $j, $offx, $offy);
#		        if ($bkgnd > $bkg[$j]) {
#		            $bkgnd = $bkg[$j] + $SKYVAL;
#		        };
                    };
		};
	    };		

     # Make sky component
            $ncomp++;
            sky ($ncomp, $bkgnd);

     # Put image fitting size and convolution box info into the galfit 
     # template file

#	    if ($ncomp == 2) {
	        $cbox = 60;
#	    } else {
#	        $cbox = int(2 * $size + 1);
#	    };

            print OUT1 ("H) $xmin $xmax $ymin $ymax # Image region to fit (xmin xmax ymin ymax) \n");
            print OUT1 ("I) $cbox   $cbox       # Size of the convolution box (x y) \n");

     # Calculate information needed to plug back into the image header
     # at the end of the fit.

	    $xmin = $xmin + $offx;  # The [xmin:xmax,ymin:ymax] of the box
	    $xmax = $xmax + $offx;  # relative to the big ACS image from which
	    $ymin = $ymin + $offy;  # the current tile was extracted.
            $ymax = $ymax + $offy;
    
            print OUT2 ("$sobj $outname-out   $offx   $offy  [$xmin:$xmax,$ymin:$ymax] $x[$i] $y[$i]\n");

            close (OUT1);

            $errno = system ("galfit $parfile");   	    # Run galfit

            if ($errno == 25600) {
		print OUT3 ("$sobj $parfile \n");
	    } else {
		print OUT4 ("$outname-out.fits \n");
	    };
#	print ("\n\n###########################\n");
#       print ("The error number is: $errno\n");
#	print ("###########################\n\n");

        };
    };


}

##############################################################################
##############################################################################

sub components {
    my $objn = @_[0]; my $centobj = @_[1], $curobj = @_[2]; 
    my $xoff = @_[3]; my $yoff = @_[4];
    my ($xpos, $ypos);
    my ($msersic, $Re_sersic, $ar_sersic, $n, $mexpdisk, $Rs);

 # Calculate the (x,y) position of the current object relative to
 # the tile in which it lives.

    $xpos = $x[$curobj] - $xoff;
    $ypos = $y[$curobj] - $yoff;

    if ($FITFUNC eq "sersic" || $centobj != $curobj) {
	$n = 1.5;
	$msersic = $mag[$curobj];
	$Re_sersic = $Re[$curobj];
	$ar_sersic = $axrat[$curobj];
    } elsif ($centobj == $curobj) {
	$n = 1.5;
	$msersic = $mag[$curobj] + 0.5;
	$ar_sersic = 0.8;
	$Re_sersic = $Re[$curobj] * 0.7;
	$mexpdisk = $mag[$curobj] + 0.5;
	$Rs = $Re[$curobj];
    }

    print OUT1 ("# Object number: $objn \n");
    print OUT1 (" 0)     sersic         #    Object type \n");
    print OUT1 (" 1) $xpos  $ypos  1 1    #   position x, y \n");
    print OUT1 (" 3) $msersic   1       #  total magnitude  \n");
    print OUT1 (" 4) $Re_sersic 1       #      R_e \n");
    print OUT1 (" 5) $n         1       #  exponent (de Vaucouleurs = 4) \n");
    print OUT1 (" 9) $ar_sersic 1       #  axis ratio (b/a) \n");
    print OUT1 ("10) $theta[$curobj] 1       #  position angle (PA) \n");
    print OUT1 (" Z) 0                  #  Output option (0 = residual, 1 = Don't subtract) \n");
    print OUT1 ("\n");

    if ($curobj == $centobj && $FITFUNC eq "BD") {
	$objn++;
        print OUT1 ("# Object number: $objn \n");
        print OUT1 (" 0)     expdisk         #    Object type \n");
        print OUT1 (" 1) $xpos  $ypos  1 1    #   position x, y \n");
        print OUT1 (" 3) $mexpdisk  1       #  total magnitude  \n");
        print OUT1 (" 4) $Rs        1       #      R_e \n");
        print OUT1 (" 9) $axrat[$curobj] 1       #  axis ratio (b/a) \n");
        print OUT1 ("10) $theta[$curobj] 1       #  position angle (PA) \n");
        print OUT1 (" Z) 0                  #  Output option (0 = residual, 1 = Don't subtract) \n");
        print OUT1 ("\n");
    };
    return ($objn);
}

sub sky {
    my $objn = @_[0]; my $bkgd = @_[1];
    
    print OUT1 ("# Object number: $objn \n");
    print OUT1 (" 0)        sky         #    Object type \n");
    print OUT1 (" 1) $bkgd     1       #  sky background \n");
    print OUT1 (" Z) 0                #  Output option (0 = residual, 1 = Don't subtract) \n");
    print OUT1 ("\n");
    print OUT1 ("================================================================================\n");

}

##############################################################################
##############################################################################


#subroutine deg2cel.precise  (version 06/09/02)
#----------
sub deg2celestial {
  #Call w/ ($ra_celestial,$dec_celestial) = deg2celestial($ra,$dec);
  #returns celestial coordinates from input of ra and dec in decimal degree
  #format.  This subroutine returns all significant digits in seconds of
  #RA and Dec; thus simpler version of deg2celestial.sub w/out roundoff.
  # NOTE: Run this subroutine in a loop, only does one transf. at a time.
  # 6/9/02: fixed bug for -1<Dec<0.

  #local variable declaration
  my $ra = @_[0];
  my $dec = @_[1];
  my $rah; my $ram; my $ras; my $hr; my $junk; my $min; my $sec;
  my $decdeg; my $decmin; my $decsec; my $sign; my $junk1; my $degree;
  my $dmin; my $dsec;
  my $ra_celestial; my $dec_celestial;

 ##RA portion
  $hr = $ra / 15.;
  ($rah,$junk)=split '\.',$hr;
  $min = 60.*("." . $junk);
  ($ram,$junk)=split '\.',$min;
  $ras = 60.*("." . $junk);

 #always want double digits
  if ($rah < 10) {$rah="0" . $rah}
  if ($ram < 10) {$ram="0" . $ram}
  if ($ras < 10) {$ras="0" . $ras}
  $ra_celestial = $rah . ":" . $ram . ":" . $ras;

 ##Dec portion
  ($degree,$junk)=split '\.',$dec;
  if ($degree >= 0.0) {
     $sign = "p";
     ($junk1,$deg)=split '\+',$degree;
  }
  if ($degree <= 0.0 && $dec < 0.0) {
     $sign = "m";
     ($junk1,$deg)=split '\-',$degree;
  }
  $dmin = 60.*("." . $junk);
  ($amin,$junk)=split '\.',$dmin;
  $asec = 60.*("." . $junk);

 #always want double digits
  if ($deg < 10) {$deg="0" . $deg}
  if ($amin < 10) {$amin="0" . $amin}
  if ($asec < 10) {$asec="0" . $asec}

 #put sign back on
  if ($sign eq "m") {$dec_celestial = "-" . $deg . ":" . $amin . ":" . $asec}
  if ($sign eq "p") {$dec_celestial = "+" . $deg . ":" . $amin . ":" . $asec}

  return ($ra_celestial,$dec_celestial);
}

#############################################################################

sub uniqname 
{

  my $ra_celest = @_[0];
  my $dec_celest = @_[1];
  my $ra1; my $ra2; my $ra3; my $ra_name;
  my $d1; my $d2; my $d3; my $dec_name;
  my $out_name; my $sign;

  ($ra1,$ra2,$ra3)=split '\:',$ra_celest;
  $ra3 = sprintf("%.2f", $ra3);      #rounds sec to 0.01
  if ($ra3 < 10.) {$ra3 = "0" . $ra3}
  $ra_name = $ra1 . $ra2 . $ra3;
  ($d1,$d2,$d3)=split '\:',$dec_celest;
  $d3 = sprintf("%.1f", $d3);        #rounds asec to 0.1
  if ($d3 < 10.) {$d3 = "0" . $d3}
  if ($d1 < 0.) {$sign="m"}
  if ($d1 >= 0.) {$sign="p"}
  $d1 = substr($d1,1,2);             #deletes sign
  $dec_name = $sign . $d1 . $d2 . $d3;
  $out_name = "GEMS" . $ra_name . $dec_name;
  return ($out_name);
}


##############################################################################

# Choose the appropriate convolution PSF.  Given a list of PSF names with
# names PSF-x-y.fits, located at position (x, y), figure out for a galaxy
# at position xgal, ygal, which PSF is the closest to this object.
#

sub psfdist
{
    my $xg = @_[0];
    my $yg = @_[1];

    my ($dum, $xp, $yp);
    my ($nearpsf, $usepsf, $psfs);
    my ($i, $dist, $neardist);

    $i = 0;
    open (PSFFILE, "psf.temp");
    while ($psfs[$i]=<PSFFILE>) {
        chomp ($psfs[$i++]);
    };

    $i = 0;
    $neardist = 1.e6;
    while ($psfs[$i]) {
        ($dum, $xp, $yp) = split '-', $psfs[$i];
        ($yp) = split '.fits', $yp;
        if (($dist = sqrt(($xp - $xg)**2 + ($yp - $yg)**2)) < $neardist) {
            $neardist = $dist; 
            $usepsf = $psfs[$i];
        }
        $i++;
    };    

    close (PSFFILE);
    return ($usepsf);
}