#!/bin/bash
#! /Users/ana/software/Ureka/Ureka/variants/common/bin/python
# SDSS DR8 data

E_DATAFILE=90
if [ $# -ge 2 ]; then image="$1"; sexfile="$2" ;else
echo "Usage: $0 image.fits .sex "
exit $E_DATAFILE; 
fi #  ASCII text file

if [ ! -e "$image" ]
then
  echo "\""$image"\" does not exist in the current dir!"
  exit $E_DATAFILE
fi
galaxy=$(basename $image);
galaxyName=$(basename $image | cut -d. -f1)
echo $galaxyName
#paramfile=$(echo $sexfile | sed 's/sex/param/g')
paramfile="mass.param"
if [ ! -e "$sexfile" ] || [ ! -e "$paramfile" ]; then
echo "\""$sexfile" or "$paramfile"\" do not exist in the current dir!"
exit $E_DATAFILE
fi

if [ -z `which sex` ]; then echo "You have to install SExtractor first."
			    exit
fi


if [ -z `which imhead` ]; then echo "You have to install WCSTools first:"
			       echo "http://tdc-www.harvard.edu/software/wcstools/index.html"
			       exit
fi
# ------------------------ S T A R T --------------------------------
mag0=$(imhead $image | grep -i magzp | awk '{print $3}')
mu0=$mag0
echo "mag0: $mag0   mu0: $mu0"
#magzp=$(echo "$exptime $pixScale $aterm $kterm $kdot $tai $reftai $airmass" | awk '{print 2.5*log($1)/log(10)+2.5*log($2*$2)/log(10)-$3-($4+$5*($6-$7))*$8}')
#magzp="24.8"
# SExtractor call to make segmentation image and determine FWHM as a crutial param for star/galaxy separation
if [ -e $galaxyName.mask.fit ]; then rm $galaxyName.mask.fit; fi
saturValue="5000"
#bkgFile="rimages.median.dat"
#if ! (test -s "$bkgFile"); then  echo "BKG keyword missing!"; exit; fi
#background=`grep -iw $galaxyName $bkgFile | awk '{print $2}' `
#background="0.003119633"
sed -e s/galaxyImage/$image/ -e s/nsigma/3/ background.py > calc_background.py
#which python
#background=$(python calc_background.py | awk '{print $2}' | awk -F= '{sum+=$2} END {print $2/NR}')
background=$(python calc_background.py | tail -1 | awk '{print $5}' | awk -F= '{print $2}')
python calc_background.py | tail -1 | awk '{print $2}' | awk -F= '{print "2:" $2}'
python calc_background.py | tail -1 | awk '{print $5}' | awk -F= '{print "5:" $2}'
echo "name:$galaxy mag0:$mag0 gain:$gain saturation:$saturation background:$background "
if [ -z "$background" ]; then echo "BKG keyword missing!"; exit; fi

# Determine background inside 512x512 mesh grid and subtract it from the image -- make galaxy.nobkg.fit image
#sed -e 's/galaxyName/'$galaxy'/g'  -e 's/magZero/'$mag0'/g' -e 's/gainValue/'$gain'/g' -e 's/saturValue/'$saturation'/g' -e 's/backgroundValue/'$background'/g' $sexfile
sed -e 's/galaxyName/'$galaxyName'/g' -e 's/magZero/'$mag0'/g' -e 's/backgroundValue/'$background'/g' -e  's/checkValue/-BACKGROUND/g' $sexfile > $galaxyName.sex
sex -c $galaxyName.sex $image # > tmp 2>&1 ## 2 Standard output for errors; 1 Standard output
#sed -e 's/galaxyName/'$galaxy'/g'  -e 's/magZero/'$mag0'/g' -e 's/gainValue/'$gain'/g' -e 's/saturValue/'$saturation'/g' nobkg$sexfile > $galaxy.sex
#sex -c $galaxy.sex $galaxy.nobkg.fit 
#rm $galaxyName.sex
# -0.0001054146
datafile="$galaxyName.cat"
if [ `cat $datafile | grep \# | wc -l` -lt 14 ]; then 
echo "ERROR: The file must be of ASCII_HEAD type in SExtractor. Minimum fields required: "
echo "       X_IMAGE, Y_IMAGE, MAG_ISO, ISOAREA_IMAGE, BACKGROUND, FLUX_RADIUS, FWHM_IMAGE"
exit
fi

#date '+%H:%M:%S'
s1=`date '+%H:%M:%S'`
comments=0
for file in `ls $datafile` 
do
  while read dataline   # Read one line at a time.
  do
     if [ ${dataline:0:1} = "#" ]; then
     tmp=`echo $dataline | awk '{print $3}'`
     let comments=$comments+1
     echo $dataline >> header
     #echo $dataline | awk '{print $2}'
     case  "$tmp" in
     "X_IMAGE")       ximage=`echo $dataline | awk '{print $2}'`;;
     "Y_IMAGE")       yimage=`echo $dataline | awk '{print $2}'`;;
     "NUMBER")        number=`echo $dataline | awk '{print $2}'`;;
     "FLAGS")         flags=`echo $dataline | awk '{print $2}'`;;
     "MAG_ISO")       magiso=`echo $dataline | awk '{print $2}'`;;
     "ISOAREA_IMAGE") isoarea=`echo $dataline | awk '{print $2}'`;;
     "BACKGROUND")    bkg=`echo $dataline | awk '{print $2}'`;;
     "FLUX_RADIUS")   fluxrad=`echo $dataline | awk '{print $2}'`;;
     "FLUX_AUTO")     fluxauto=`echo $dataline | awk '{print $2}'`;;
     "MAG_AUTO")      magauto=`echo $dataline | awk '{print $2}'`;;
     "FWHM_IMAGE")    fwhm=`echo $dataline | awk '{print $2}'`;;
     "X_WORLD")       xworld=`echo $dataline | awk '{print $2}'`;;
     "Y_WORLD")       yworld=`echo $dataline | awk '{print $2}'`;;
     "XMIN_IMAGE")    xmin=`echo $dataline | awk '{print $2}'`;;
     "XMAX_IMAGE")    xmax=`echo $dataline | awk '{print $2}'`;;
     "YMIN_IMAGE")    ymin=`echo $dataline | awk '{print $2}'`;;
     "YMAX_IMAGE")    ymax=`echo $dataline | awk '{print $2}'`;;
     "A_IMAGE")       Aimg=`echo $dataline | awk '{print $2}'`;;
     "B_IMAGE")       Bimg=`echo $dataline | awk '{print $2}'`;;
     "THETA_IMAGE")   THETAimg=`echo $dataline | awk '{print $2}'`;; 
     * );;
     esac
     else break
     fi
  done < $file 
done 
#}   
  if [ -z "$ximage" ] || [ -z "$yimage" ] || [ -z "$magiso" ] || [ -z "$isoarea" ] || [ -z "$bkg" ] || [ -z "$fluxrad" ] || [ -z "$fluxauto" ] || [ -z "$magauto" ]  || [ -z "$fwhm" ] || [ -z "$xworld" ] || [ -z "$yworld" ] || [ -z "$xmin" ] || [ -z "$xmax" ] || [ -z "$ymin" ] || [ -z "$ymax" ] || [ -z "$Aimg" ] || [ -z "$Bimg" ] || [ -z "$THETAimg" ] || [ -z "$fwhm" ]; then
  echo "You are missing one of the obligated parameters - check the file again:";
  echo "X_IMAGE, Y_IMAGE, MAG_ISO, ISOAREA_IMAGE, BACKGROUND, FLUX_RADIUS, FLUX_AUTO, MAG_AUTO, FWHM_IMAGE, X_WORLD, Y_WORLD, XMIN_IMAGE, XMAX_IMAGE, YMIN_IMAGE, YMAX_IMAGE, A_IMAGE, B_IMAGE, THETA_IMAGE, FWHM_IMAGE"
  exit
  fi
# ------------------------- FIND a galaxy in the sex image according to coordinates provided by SDSS ----------------------------
# 2MASS pixel scale: 1.
radius="0.00083" # 3 [pix] * 1. ["] = {seeing ["]} = seeing["]/3600. =seeing [deg]
alpha=$(imhead $image | grep -w CRVAL1 | awk '{printf "%.8f\n", $3}')
delta=$(imhead $image | grep -w CRVAL2 | awk '{printf "%.8f\n", $3}')
if [ -z "$alpha" ] || [ -z "$delta" ]; then echo "Image name MUST contain the name of the galaxy that is to be found in *$photfile* file";
			     exit;
fi
echo "(ra:dec)=($alpha:$delta)"

imageCoords=""
while [ -z "$imageCoords" ] && [ `echo "$imageCoords" | wc -l` -le 1 ]; do
# Search within seeing radius=0.00033 if there is an entry (galaxy) in the SExtractor .cat file: ALPHA=X_WORLD+/-SEEING and DELTA=Y_WORLD+/-SEEING
imageCoords=`awk -v a="$alpha" -v d="$delta" -v r="$radius" -v c="$comments" -v xw=$xworld -v yw=$yworld -v x=$ximage -v y=$yimage -v xl=$xmin -v xu=$xmax -v yl=$ymin -v yu=$ymax -v aimg=$Aimg -v bimg=$Bimg -v t=$THETAimg -v m=$magiso -v bkg=$bkg -v fauto=$fluxauto -v mauto=$magauto  -v num=$number -v fwhm=$fwhm '{
	     if (NR > c && $xw > a-r && $xw < a+r && $yw > d-r && $yw < d+r) printf "%.3f   %.3f   %d  %d  %d  %d  %.3f   %.3f   %.1f   %.3f   %.3f   %.3f   %.3f   %5d %.3f\n", $x, $y, $xl, $xu, $yl, $yu, $aimg, $bimg, $t, $m, $bkg, $fauto, $mauto, $num, $fwhm}' $datafile`
# If not, check if X_WORLD and Y_WORLD coordinates are switched X_WORLD <-> YWORLD incidently 
if [ -z "$imageCoords" ]; then 
imageCoords=`awk -v a="$alpha" -v d="$delta" -v r="$radius" -v c="$comments" -v yw=$xworld -v xw=$yworld -v x=$ximage -v y=$yimage -v xl=$xmin -v xu=$xmax -v yl=$ymin -v yu=$ymax -v aimg=$Aimg -v bimg=$Bimg -v t=$THETAimg -v m=$magiso -v bkg=$bkg -v fauto=$fluxauto -v mauto=$magauto  -v num=$number -v fwhm=$fwhm '{
	     if (NR > c && $xw > a-r && $xw < a+r && $yw > d-r && $yw < d+r) printf "%.3f   %.3f   %d  %d  %d  %d  %.3f   %.3f   %.1f   %.3f  %.3f   %.3f   %.3f   %5d %.3f\n", $x, $y,  $xl, $xu, $yl, $yu, $aimg, $bimg, $t, $m, $bkg, $fauto, $mauto, $num, $fwhm}' $datafile`
fi
# If nothing is found both for (ALPHA,DELTA)=(X_WORLD,Y_WORLD)+/-SEEING and (ALPHA,DELTA)=(Y_WORLD,X_WORLD)+/-SEEING duplicate searching radius until a match is found
radius=$(echo "scale=6; $radius + $radius" | bc)
done
if [ -z "$imageCoords" ] || [ `echo "$imageCoords" | wc -l` -ne 1 ]; then 
echo "A unique match between equ. (J2000) and SExtractor world coords NOT FOUND!"; echo "Check both photo and .cat files."; exit
fi     

# --------- FIND a rectangle that encompasses most of the galaxy flux -------------
xc=`echo $imageCoords | awk '{print $1}'`
xcenter=$xc 
yc=`echo $imageCoords | awk '{print $2}'`
ycenter=$yc
num=`echo $imageCoords | awk '{print $(NF-1)}'`
pixScale="1"
fwhm=`echo $imageCoords | awk -v pixScale=$pixScale '{print $NF*pixScale}'`
echo "Object #: xc[pix]  yc[pix]     FWHM[\"]"
echo "       "$num":" $xcenter $ycenter "   " $fwhm
#echo $imageCoords
echo "       ISO               AUTO                APER"
cat $datafile | grep " $xc " | grep " $yc " | awk '{print $8"+/-"$9,$10"+/-"$11,$12"+/-"$13}'
xlow=`echo $imageCoords | awk '{print $3}'`
xup=`echo $imageCoords | awk '{print $4}'`
ylow=`echo $imageCoords | awk '{print $5}'`
yup=`echo $imageCoords | awk '{print $6}'`
A=`echo $imageCoords | awk '{print $7}'`
B=`echo $imageCoords | awk '{print $8}'`
THETA=`echo $imageCoords | awk '{print $9}'`
magnitude=`echo $imageCoords | awk '{print $10}'`
BKGimg=`echo $imageCoords | awk '{print $11}'`
#echo "SExtractor background: $BKGimg"
flux=`echo $imageCoords | awk '{print $12}'`
mag=`echo $imageCoords | awk '{print $13}'`
nn=`echo $imageCoords | awk '{print $14}'` # The serial number of the galaxy in the image
ellip=`echo "scale=4; 1 - $B/$A" | bc -l`
##if [ $(echo ${posAngle:0:1}) = "-" ]; then 
# P.A. goes from -90 to 90 
##posAngle=`echo "scale=3; $THETA + 90." | bc`; 
# echo "290" | awk '{printf "%d\n", $1/90}'
##else posAngle=$THETA; fi

#posAngle=`echo "scale=3; $THETA + 90." | bc`; # PROMENA
posAngle="$THETA"
testAngle=`echo "$THETA/1" | bc`;
if [ "$testAngle" -gt 90 ] && [ "$testAngle" -lt 270 ]; then posAngle=`echo "$THETA-180" | bc -l`; fi
if [ "$testAngle" -gt 270 ]; then posAngle=`echo "$THETA-360" | bc -l`; fi
#posAngle=$(echo $THETA | awk 'function abs(t){return (((t < 0.0) ? -t : t) + 0.0)}{print 90-abs($1)}')
BA=$(echo "scale=1; ($B/$A)*100." | bc | cut -d'.' -f1)
xsize=$(imhead $image | grep NAXIS1 | awk '{print $3}')
ysize=$(imhead $image | grep NAXIS2 | awk '{print $3}')
if [ -z "$xsize" ]; then 
xsize=`imsize $image | cut -d'x' -f3`
ysize=`imsize $image | cut -d'x' -f4 | cut -d'p' -f1`
fi 
diffxl=$(echo "scale=0; $xc/1.0-$xlow/1.0" | bc)
diffxu=$(echo "scale=0; $xup/1.0-$xc/1.0" | bc)
if [ `echo $[$diffxu-$diffxl]` -lt 2 ] || [ `echo $[$diffxu-$diffxl]` -gt 2 ]; then 
   if [ $diffxl -gt $diffxu ]; then xlow=$(echo "scale=0; $xc/1.0+$xc/1.0-$xup/1.0"| bc); else xup=$(echo "scale=0; $xc/1.0+$xc/1.0-$xlow/1.0"| bc); fi
fi   
diffyl=$(echo "scale=0; $yc/1.0-$ylow/1.0" | bc)
diffyu=$(echo "scale=0; $yup/1.0-$yc/1.0" | bc)
if [ $(echo $[$diffyu-$diffyl]) -lt 2 ] || [ $(echo $[$diffyu-$diffyl]) -gt 2 ]; then 
   if [  $diffyl -gt $diffyu ]; then ylow=$(echo "scale=0; $yc/1.0+$yc/1.0-$yup/1.0"| bc); else yup=$(echo "scale=0; $yc/1.0+$yc/1.0-$ylow/1.0"| bc); fi
fi   
limit=250
if [ $xlow -le $limit ]; then xl=1; xu=$[$xup+$xlow]; else 
		      if [ `echo $[$xup+$limit]` -gt $xsize ]; then xu=$xsize; xl=$[$xlow - ($xsize-$xup)]; 
		      else let xu=$xup+$limit; let xl=$xlow-$limit;
		      fi
fi
if [ $ylow -le $limit ]; then yl=1; let yu=$[$yup+$ylow]; else 
		      if [ `expr $yup+$limit | bc` -gt $ysize ]; then yu=$ysize; yl=$[$ylow - ($ysize-$yup)]; 
		      else let yu=$yup+$limit; let yl=$ylow-$limit;
		      fi
fi
# Changing xc and yc to correspond to the cut off image
xc=`echo "scale=3; ($xu - $xl)/2." | bc`
yc=`echo "scale=3; ($yu - $yl)/2." | bc`

# starting NEW as of 25 July 2014 to spped up things
#skymag=`echo "scale=6; $magzp-$extinction*$airmass" | bc | awk -v exptime=$exptime '{print $1+2.5*log(exptime)/log(10)+2.5*log(0.16)/log(10)}'` #mag0
# skymag == surface brightness zero point adding 2.5log10((")^2)=2.5log10(pixScale^2)  pixScale=0.396"
skymag=$mag0
skyflux="1"  # refer
bias="0" # zerolevel
echo "... Running ELLIPSE :: refer = $skyflux    mu0 = $skymag    zerolevel = $bias"
echo "major   minor axis     eccentricity   PA"
echo "$A         $B           $ellip        $posAngle"

#background=`imhead $smallImage | grep -iw bkg | awk '{print $3}'` # Sextractor estimation
# imcopy r.fits[0][1:1854,1:1122] NGC3198.short.fit
sed -e 's/galaxyName/'$galaxyName'/g' -e 's/galaxy/'$galaxy'/g' -e 's/skyMag/'$skymag'/g' -e 's/ecc/'$ellip'/g' -e 's/theta/'$posAngle'/g' -e 's/majoraxis/'$A'/g' -e 's/minoraxis/'$B'/g' -e 's/xc/'$xcenter'/g' -e 's/yc/'$ycenter'/g' -e 's/num/'$num'/g' -e 's/skyFlux/'$skyflux'/g' -e 's/bias/'$bias'/g' dr7_isophote.py > doEllipse_$galaxyName.py # PROMENA
#/Users/ana/software/Ureka/Ureka/variants/common/bin/python doEllipse_$galaxyName.py
exit
python doEllipse_$galaxyName.py
# ending NEW
exit # EXIT
# ------------------------------------ FINDING STARS for IRAF --------------------------------------------------
     awk -v cf="$fwhm" -v c="$comments" '{if (NR > c && $cf > 0) print $cf}' $datafile > tmp
     stdevFWHM=`stdev tmp`
     medianFWHM=`median tmp`
#     echo $stdevFWHM $medianFWHM
     awk -v cm="$magiso" -v c="$comments" '{if (NR > c && $cm > 0) print $cm}' $datafile > tmp
     stdevMAG=`stdev tmp` # Popravi stdev!!!
     averageMAG=`average tmp`
#     echo $stdevMAG $averageMAG
     maglim=$(echo "scale=6; $averageMAG - $stdevMAG/4. ;" | bc)
     fwhmlim=$(echo "scale=6; $medianFWHM + 1.;" | bc)

cp header tmp1
awk -v cx=$ximage -v cy=$yimage -v cm=$magiso -v fwhm=$fwhmlim -v f=$flags -v frad=$fluxrad -v cfwhm=$fwhm  -v c=$comments -v mag=$maglim '{
 if (NR > c && $f == 0 && $cm < mag && fwhm > 2*$frad && $cfwhm < fwhm && $cfwhm > $frad) print $0}' $datafile >> tmp1
if [ `cat tmp1 | grep -v "\#" | wc -l` -lt 5 ]; then echo "1: Not enough stars found"; exit; fi
# ----------------- CHECKING if there is less than 5 PSF stars and if so, reduce the radius from 30 ---------------
psfrad=35
if [ -e galaxyCoords ]; then rm galaxyCoords; fi
cp header galaxyCoords
while [ `cat galaxyCoords | grep -v "\#" | wc -l` -lt 5 ] && [ $psfrad -ge 10 ]; do
if [ -e galaxyCoords ]; then rm galaxyCoords; fi
let psfrad=$psfrad-5	
 for file in `ls tmp1` 
 do
  while read dataline   # Read one line at a time.
  do
  if [ ${dataline:0:1} != "#" ]; then
  let psf=$psfrad+5
  xcoord=$(echo $dataline | awk -v cx=$ximage '{print $cx}')
  ycoord=$(echo $dataline | awk -v cy=$yimage '{print $cy}') 
  index=`awk -v x=$xcoord -v y=$ycoord -v cx=$ximage -v cy=$yimage -v c=$comments -v xs=$xsize -v ys=$ysize -v psf=$psfrad  -v psfout=$psf 'function abs(t){return (((t < 0.0) ? -t : t) + 0.0)} {
  if (NR > c && (abs(x-$cx) < psfout && abs(y-$cy) < psfout) && (x-$cx != 0 || y-$cy !=0) ) print x,"\t",y}' $datafile | wc -l`
#  if (NR > c && (abs(x-$cx) < psf+5 && abs(y-$cy) < psf+5) && (x-$cx != 0 || y-$cy !=0) && (x-psf>0 && x+psf<xs && y-psf>0 && y+psf<ys)) print x,"\t",y}' $datafile | wc -l`
  if [ $index -eq 0 ]; then 
  fits=`awk -v x=$xcoord -v y=$ycoord -v cx=$ximage -v cy=$yimage -v c=$comments -v xs=$xsize -v ys=$ysize -v psf=$psf -v psfout=$psf 'function abs(t){return (((t < 0.0) ? -t : t) + 0.0)} {
  if (NR > c && (abs(x-$cx) > psfout && abs(y-$cy) > psfout && x-psf>0 && x+psf<xs && y-psf>0 && y+psf<ys) && (x-$cx != 0 || y-$cy !=0) ) print x,"\t",y}' $datafile | head -1`
  if [ ! -z "$fits" ]; then cat tmp1 | grep -w "$xcoord" | grep -w "$ycoord" >> galaxyCoords; fi
  fi 
  else echo $dataline >> galaxyCoords	
  fi
  done < $file
 done
done
if [ `cat galaxyCoords | grep -v "\#" | wc -l` -lt 5 ]; then echo "Not enough stars found"; exit; fi

# --------------------------------------- PREPARING and RUNNING iraf script ------------------------------------------
echo "... Running DAOPHOT package tasks.."
if [ -e $image.mag.0 ]; then rm $image.*; fi
# Edit image header (gain, rdnoise, objname, magzp, bkg): iraf.hedit 
# PSF fitting with moffat15 function: iraf.daophot.digiphot.[pstselect,psf,seepsf]
#sed -e 's/galaxyName/'$galaxyName'/g' -e 's/galaxyFile/'$image'/g' -e 's/magZero/'$magzp'/g' -e 's/galaxyPsfRad/'$psfrad'/g' phot.py > doPhotometry.py 
  if [ -e $galaxyName.psf.fit ]; then rm $galaxyName.psf.fit; fi 
if [ -e $galaxyName.short.fit ]; then rm $galaxyName.short.fit; fi 
if [ -e $galaxyName.mask.fit ]; then rm $galaxyName.mask.fit; fi
#if [ -e $galaxyName.cut.fit ]; then rm $galaxyName.cut.fit; fi
sed -e 's/galaxyName/'$galaxyName'/g' -e 's/galaxyFile/'$image'/g' -e 's/gainValue/'$gain'/g' -e 's/rdnoiseValue/'$rdnoise'/g' -e 's/amassValue/'$airmass'/g' -e 's/magZero/'$magzp'/g' -e 's/background/'$BKGimg'/g' -e 's/galaxyPsfRad/'$psfrad'/g' -e 's/xlow/'$xl'/g' -e 's/xup/'$xu'/g' -e 's/ylow/'$yl'/g' -e 's/yup/'$yu'/g' photpsfcutNOBKG.py > doAll.py
#sed -e 's/galaxyName/'$galaxyName'/g' -e 's/galaxyFile/'$image'/g' -e 's/gainValue/'$gain'/g' -e 's/rdnoiseValue/'$rdnoise'/g' -e 's/amassValue/'$airmass'/g' -e 's/magZero/'$magzp'/g' -e 's/background/'$BKGimg'/g' -e 's/galaxyPsfRad/'$psfrad'/g' -e 's/xlow/'$xl'/g' -e 's/xup/'$xu'/g' -e 's/ylow/'$yl'/g' -e 's/yup/'$yu'/g' photpsfcut.py > doAll.py
python doAll.py
# ----------------------------- Mask image making -------------------------------------
# Mask image $galaxy.pl.fit - all pixels corresp to galaxy and background are set to 0.
if [ -e allzero.fit ]; then rm allzero.fit; fi
#if [ -e $galaxyName.large.mask.fit ]; then rm $galaxyName.large.mask.fit ; fi  # PROMENA
sed -e 's/galaxyName/'$galaxyName'/g' -e 's/galnumber/'$nn'/g' final.py > $galaxy.final.py # PROMENA
python $galaxy.final.py; rm $galaxy.final.py  
mv allzero.fit $galaxyName.pl.fit; 

if [ "$dim" -eq "1" ] || [ "$dim" -eq "3" ]; then 
########################################################################################################################
#------------------------------------- ELLIPSE fitting procedure -----------------------------------------
#  Making final galaxy image - all other objects masked 
# 1) expr (galaxyName.mask.fit >= 1 && galaxyName.mask.fit != galaxy_number) ? 0 : 1 -> zero.fit 
# 2) zero.fit * galaxyName.cut.fit = allzero.fit (all objects exept galaxy and background set to zero)
#echo "... Preparing galaxy image for ELLIPSE task."
# ----------------------------- Background subtraction -------------------------------------
# Image $galaxyName.cut.fit without background is created
#background=`imhead $galaxy.short.fit | grep -iw sky | awk '{print 1000.+$3}'` # SDSS header information, 1000 is the 'soft-bias'
#if [ -z "$background" ]; then background=$BKGimg; fi

#awk -v b=$bkg '!($0 ~ /#/ ){print $b}' $galaxyName.cat > back
#background=`median back`
#if [ -e $galaxy.cut.fit ]; then rm $galaxy.cut.fit; fi
#sed -e 's/galaxyName/'$galaxy'/g' -e 's/const/'$background'/g' imarith.py > $galaxy.back.py  # PROMENA!!!
#python $galaxy.back.py; rm $galaxy.back.py 
 cp $galaxy.short.fit $galaxy.cut.fit # When image without bkg is needed
# mag = mag0 - 2.5 log[(intens - zerolevel)/refer]; mag0 = -aa - kk*airmass; refer = background (DN); zerolevel=(softbias)= 1000 (DN) za SDSS
# !!! U photpsfcut.py sam dodala da isece manju sliku od "$galaxyName.nobkg.fit" = slika bez pozadine i to je photpsfcutNOBKG.py
skymag=`echo "scale=6; $magzp-$extinction*$airmass" | bc | awk -v exptime=$exptime '{print $1+2.5*log(exptime)/log(10)+2.5*log(0.16)/log(10)}'` #mag0 
skyflux="1"  # refer
bias="0" # zerolevel
echo "... Running ELLIPSE :: refer = $skyflux    mag0 = $skymag    zerolevel = $bias"
nobkgImage="$galaxyName.cut.fit"
#background=`imhead $smallImage | grep -iw bkg | awk '{print $3}'` # Sextractor estimation
sed -e 's/galaxyName/'$galaxyName'/g' -e 's/galaxyNew/'$nobkgImage'/g' -e 's/skyMag/'$skymag'/g' -e 's/ecc/'$ellip'/g' -e 's/theta/'$posAngle'/g' -e 's/majoraxis/'$A'/g' -e 's/minoraxis/'$B'/g' -e 's/xc/'$xc'/g' -e 's/yc/'$yc'/g' -e 's/skyFlux/'$skyflux'/g' -e 's/bias/'$bias'/g' isophote.py > doEllipse_$galaxyName.py # PROMENA
# The change : use the whole image instead the cut one 
#sed -e 's/galaxyName/'$galaxyName'/g' -e 's/galaxyFile/'$smallImage'/g' -e 's/galaxyNew/'$nobkgImage'/g' -e 's/skyMag/'$skymag'/g' -e 's/ecc/'$ellip'/g' -e 's/theta/'$posAngle'/g' -e 's/majoraxis/'$A'/g' -e 's/minoraxis/'$B'/g' -e 's/xc/'$xcenter'/g' -e 's/yc/'$ycenter'/g' -e 's/skyFlux/'$skyflux'/g' -e 's/bias/'$bias'/g' -e 's/bkg/'$background'/g' isophote.py > doEllipse_$galaxyName.py # PROMENA
if [ -e $galaxyName.tab ]; then rm $galaxyName.tab; fi
if [ -e mask.fit ]; then rm mask.fit; fi
if [ -e noback.fit ]; then rm noback.fit; fi
python doEllipse_$galaxyName.py > $galaxyName.txt

if [ -e temporary2 ]; then rm temporary*; fi
oldf=1;i=0;c=1;j=0;oldline="#"
while read line
do  
 let i=$i+1
 case $(echo ${line:0:1}) in 
  '#') f=1 ;;
  *) f=0 ;;
 esac
if [ "$oldf" -eq 0 ] && [ "$f" -eq 1 ]; then 
    if [ $i -gt 15 ]; then 
    let j=$j+1
    awk -v low=$(echo "$c-1" | bc) -v high=$i '(NR > low && NR < high && high > 1){ sub(/\#/, ""); print}' $galaxyName.txt > temporary$j
    c=$i 
    fi
 fi
oldf="$f"
done < $galaxyName.txt 
paste temporary[2-$j] > $galaxyName.dat
rm temporary* $galaxyName.txt $galaxyName.tab 
sed -i 's/\t\t*$//;/^$/d' $galaxyName.dat
fi 
if [ "$dim" -eq "2" ] || [ "$dim" -eq "3" ]; then 
###################################################################################################################
#-------------------------------------------  GALFIT fitting procedure -----------------------------------------------
if [ -e "galfit.01" ] || [ -e "galfit.02" ]; then rm galfit.0*; fi
xgal=`echo $[$xu-$xl]`; ygal=`echo $[$yu-$yl]`; 
sed -e 's/galaxyName/'$galaxyName'/g' -e 's/xsize/'$xgal'/g' -e 's/ysize/'$ygal'/g' -e 's/magZero/'$magzp'/g' -e 's/xcenter/'$(echo "scale=3; $xgal/2." | bc)'/g' -e 's/ycenter/'$(echo "scale=3; $ygal/2." | bc)'/g' -e 's/majorAxis/'$A'/g' -e 's/minorAxis/'$(echo "scale=3; $B/4." | bc)'/g' -e 's/minor2major/'$(echo "scale=3; $B/$A" | bc)'/g' -e 's/posAngle/'$(echo "scale=3; $THETA + 90." | bc)'/g' -e 's/magIso/'$magnitude'/g' -e 's/bkgValue/'$BKGimg'/g' galfit.ser.feedme > $galaxyName.ser.feedme
sed -e 's/galaxyName/'$galaxyName'/g' -e 's/xsize/'$xgal'/g' -e 's/ysize/'$ygal'/g' -e 's/magZero/'$magzp'/g' -e 's/xcenter/'$(echo "scale=3; $xgal/2." | bc)'/g' -e 's/ycenter/'$(echo "scale=3; $ygal/2." | bc)'/g' -e 's/majorAxis/'$A'/g' -e 's/minorAxis/'$(echo "scale=3; $B/4." | bc)'/g' -e 's/minor2major/'$(echo "scale=3; $B/$A" | bc)'/g' -e 's/posAngle/'$(echo "scale=3; $THETA + 90." | bc)'/g' -e 's/magIso/'$magnitude'/g' -e 's/bkgValue/'$BKGimg'/g' galfit.dev.feedme > $galaxyName.dev.feedme
sed -e 's/galaxyName/'$galaxyName'/g' -e 's/xsize/'$xgal'/g' -e 's/ysize/'$ygal'/g' -e 's/magZero/'$magzp'/g' -e 's/xcenter/'$(echo "scale=3; $xgal/2." | bc)'/g' -e 's/ycenter/'$(echo "scale=3; $ygal/2." | bc)'/g' -e 's/majorAxis/'$A'/g' -e 's/minorAxis/'$(echo "scale=3; $B/4." | bc)'/g' -e 's/minor2major/'$(echo "scale=3; $B/$A" | bc)'/g' -e 's/posAngle/'$(echo "scale=3; $THETA + 90." | bc)'/g' -e 's/magIso/'$magnitude'/g' -e 's/bkgValue/'$BKGimg'/g' galfit.exp.feedme > $galaxyName.exp.feedme
sed -e 's/galaxyName/'$galaxyName'/g' -e 's/xsize/'$xgal'/g' -e 's/ysize/'$ygal'/g' -e 's/magZero/'$magzp'/g' -e 's/xcenter/'$(echo "scale=3; $xgal/2." | bc)'/g' -e 's/ycenter/'$(echo "scale=3; $ygal/2." | bc)'/g' -e 's/majorAxis/'$A'/g' -e 's/minorAxis/'$(echo "scale=3; $B/4." | bc)'/g' -e 's/minor2major/'$(echo "scale=3; $B/$A" | bc)'/g' -e 's/posAngle/'$(echo "scale=3; $THETA + 90." | bc)'/g' -e 's/magIso/'$magnitude'/g' -e 's/bkgValue/'$BKGimg'/g' galfit.serexp.feedme > $galaxyName.serexp.feedme
sed -e 's/galaxyName/'$galaxyName'/g' -e 's/xsize/'$xgal'/g' -e 's/ysize/'$ygal'/g' -e 's/magZero/'$magzp'/g' -e 's/xcenter/'$(echo "scale=3; $xgal/2." | bc)'/g' -e 's/ycenter/'$(echo "scale=3; $ygal/2." | bc)'/g' -e 's/majorAxis/'$A'/g' -e 's/minorAxis/'$(echo "scale=3; $B/4." | bc)'/g' -e 's/minor2major/'$(echo "scale=3; $B/$A" | bc)'/g' -e 's/posAngle/'$(echo "scale=3; $THETA + 90." | bc)'/g' -e 's/magIso/'$magnitude'/g' -e 's/bkgValue/'$BKGimg'/g' galfit.devexp.feedme > $galaxyName.devexp.feedme
if [ -e $image.mag.1 ]; then rm $image.* $galaxyName.sex; fi

galfit $galaxyName.serexp.feedme; if [ ! -e $galaxyName.serexp.fits ] || [ ! -z "$(imhead $galaxyName.serexp.fits[2] | grep \*)" ]; then # Fit didn't crash,and is meaningfull 
				     echo "Trying DEV+EXP profile ..."; galfit $galaxyName.devexp.feedme; 
				       if [ -e $galaxyName.devexp.fits ] && [ -z "$(imhead $galaxyName.devexp.fits[2] | grep \*)" ]; then # Fit successful and there are no suspicious output params    
				       mv $galaxyName.devexp.fits $galaxyName.res.fits; mv galfit.01 $galaxyName.galfit; # Everything OK with deVaucouleurs+Exponential fit
				       else # Check which component causes the problem
				          if [ -e $galaxyName.devexp.fits ]; then # 
				            echo "DEV+EXP::There are some suspicious output params"; 
					    if [ ! -z "$(imhead $galaxyName.devexp.fits[2] | grep 1_ | grep \*)" ]; then # Sth wrong with deVaucouleurs component
					    echo "Trying pure Exponential profile ..."
					    rm galfit.01; galfit $galaxyName.exp.feedme; if [ -e $galaxyName.exp.fits ]; then mv $galaxyName.exp.fits $galaxyName.res.fits; mv galfit.01 $galaxyName.galfit; fi
					    fi
                                               if [ ! -z "$(imhead $galaxyName.devexp.fits[2] | grep 2_ | grep \*)" ]; then # Sth wrong with Exponential component
					       echo "Trying pure Sersic profile ..."					
					       rm galfit.01; galfit $galaxyName.ser.feedme; if [ -e $galaxyName.ser.fits ]; then mv $galaxyName.ser.fits $galaxyName.res.fits; mv galfit.01 $galaxyName.galfit; fi
					       fi 
					  else  echo "Trying pure Sersic profile ..."; # Both SER+EXP and DEV+EXP failed  
					        galfit $galaxyName.ser.feedme; if [ -e $galaxyName.ser.fits ]; then 
												mv $galaxyName.ser.fits $galaxyName.res.fits; mv galfit.01 $galaxyName.galfit; 
					      							else galfit $galaxyName.dev.feedme; if [ -e $galaxyName.dev.fits ]; then  mv $galaxyName.dev.fits $galaxyName.res.fits; mv galfit.01 $galaxyName.galfit; fi
												galfit $galaxyName.exp.feedme; if [ -e $galaxyName.exp.fits ]; then 
																chi1=`imhead $galaxyName.res.fits[2] | grep -i chi2nu | awk '{print $3}' | cut -d. -f1`; 
																chi2=`imhead $galaxyName.exp.fits[2] | grep -i chi2nu | awk '{print $3}' | cut -d. -f1`;
																    if [ $chi1 -eq "0" ]; then if [ $chi2 -eq "0" ]; then chi1=`imhead $galaxyName.res.fits[2] | grep -i chi2nu | awk '{print $3}' | cut -d. -f2`; 
																                           				  chi2=`imhead $galaxyName.exp.fits[2] | grep -i chi2nu | awk '{print $3}' | cut -d. -f2`; 	   
																						          if [ $chi1 -gt $chi2 ]; then mv $galaxyName.exp.fits $galaxyName.res.fits; mv galfit.01 $galaxyName.galfit; fi   
															                                       fi	  
																    else if [ $chi2 -eq "0" ]; then mv $galaxyName.exp.fits $galaxyName.res.fits; mv galfit.01 $galaxyName.galfit;  
																         else if [ $chi1 -gt $chi2 ]; then mv $galaxyName.exp.fits $galaxyName.res.fits; mv galfit.01 $galaxyName.galfit; fi
																	 fi
																    fi	 
																fi												
									       fi		
					  fi							
				       fi        
				  else mv $galaxyName.serexp.fits $galaxyName.res.fits; mv galfit.01 $galaxyName.galfit;  # No suspicious output params - accept the fit	
				  fi	
rm $galaxyName.s*.feedme $galaxyName.d*.feedme $galaxyName.e*.feedme $galaxyName.galfit				  
rm $galaxyName.galfit
if [ -e "galfit.01" ] || [ -e "galfit.02" ]; then rm galfit.0*; fi
if [ ! -e $galaxyName.res.fits ]; then echo "Check initial parameters and image $galaxyName.short.fit"; fi
# ------------------------------------------------------ END ------------------------------------------------------------
fi ######################################################################################################################
rm $galaxyName*mask.fit $galaxyName.short.fit $galaxyName.nobkg.fit  # $galaxyName.cut.fit
if [ -e tmp1 ]; then rm tmp1; fi
if [ -e tmp ]; then rm tmp; fi
if [ -e header ]; then rm header; fi
if [ -e $image.mag.0 ] || [ -e $image.mag.1 ]; then rm $image.mag.* ; fi
if [ -e $image.pst.1 ] || [ -e $image.pst.2 ]; then rm $image.pst.* ; fi
if [ -e $image.psg.1 ]; then rm $image.psg.* ; fi
if [ -e $image.psf.1.fits ]; then rm $image.psf.*.fits ; fi
#if [ $s2 -gt $s1 ]; then time=$[$(echo "scale=1; $s2" | bc)-$(echo "scale=1; $s1" | bc)]; else time=$[60-$(echo "scale=1; $s1" | bc)+$(echo "scale=1; $s2" | bc)]; fi
#echo "Input: A=$A  B=$B  B/A=$BA%  THETA=$THETA  MAG=$magnitude  ELLIP=$ellip"
#echo "Elapsed time: $time s"
echo "GALAXY :: $galaxyName"
echo "START: $s1"
echo " END : `date '+%H:%M:%S'`"

