#!/bin/sh

# Construct a star stamp to build a Point Spread Function (PSF). This
# script will consider a center position and then it will crop the original
# image around that center with the specified size. After that, it will
# compute the radial profile to obtain a normalization value. Finally, the
# output stamp will be normalized by dividing by the normalization value.
# More options like providing a segmentation map are also available. The
# main goal of this script is to generate stamps to combine them and get an
# extended PSF from the stamp of stars. Run with `--help', or see
# description under `print_help' (below) for more.
#
# Original author:
#   Raul Infante-Sainz <infantesainz@gmail.com>
# Contributing author(s):
#   Mohammad Akhlaghi <mohammad@akhlaghi.org>
#   Carlos Morales-Socorro <cmorsoc@gmail.com>
#   Sepideh Eskandarlou <sepideh.eskandarlou@gmail.com>
# Copyright (C) 2021-2022 Free Software Foundation, Inc.
#
# Gnuastro is free software: you can redistribute it and/or modify it under
# the terms of the GNU General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option)
# any later version.
#
# Gnuastro is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for
# more details.
#
# You should have received a copy of the GNU General Public License along
# with Gnuastro. If not, see <http://www.gnu.org/licenses/>.


# Exit the script in the case of failure
set -e

# To avoid floating points like '23,45' instead of '23.45'.
export LANG=C





# Default option values (can be changed with options on the command-line).
hdu=1
quiet=""
center=""
keeptmp=0
output=""
tmpdir=""
segment=""
axisratio=1
normradii=""
sigmaclip=""
stampwidth=""
normop="median"
positionangle=0
corewidth="3,3"
version=0.17
scriptname=astscript-psf-stamp





# Output of `--usage' and `--help':
print_usage() {
    cat <<EOF
$scriptname: run with '--help' for list of options
EOF
}

print_help() {
    cat <<EOF
Usage: $scriptname [OPTION] FITS-files

This script is part of GNU Astronomy Utilities $version.

This script will consider a center position to make a normalized stamp. It
is intendeed to generate several stamps and combine them to obtain a PSF.

For more information, please run any of the following commands. In
particular the first contains a very comprehensive explanation of this
script's invocation: expected input(s), output(s), and a full description
of all the options.

     Inputs/Outputs and options:           $ info $scriptname
     Full Gnuastro manual/book:            $ info gnuastro

If you couldn't find your answer in the manual, you can get direct help from
experienced Gnuastro users and developers. For more information, please run:

     $ info help-gnuastro

$scriptname options:
 Input:
  -h, --hdu=STR           HDU/extension of all input FITS files.
  -O, --mode=STR          Coordinates mode ('wcs' or 'img').
  -c, --center=FLT,FLT    Center coordinates of the object.
  -W, --stampwidth=INT    Width of the stamp in pixels.
  -n, --normradii=INT,INT Minimum and maximum radii (in pixels)
                          for computing the normalization value.
  -S, --segment=STR       Output of Segment (with OBJECTS and CLUMPS).
  -w, --corewidth=INT     Area width of the central object in pixels for unmasking.
  -N, --normop=STR        Operator for computing the normalization value
                          (mean, sigclip-mean, etc.).
  -Q, --axisratio=FLT     Axis ratio for ellipse maskprofile (A/B).
  -p, --positionangle=FLT Position angle for ellipse mask profile.
  -s, --sigmaclip=FLT,FLT Sigma-clip multiple and tolerance.

 Output:
  -o, --output            Output table with the radial profile.
  -t, --tmpdir            Directory to keep temporary files.
  -k, --keeptmp           Keep temporal/auxiliar files.

 Operating mode:
  -?, --help              Print this help list.
      --cite              BibTeX citation for this program.
  -q, --quiet             Don't print any extra information in stdout.
  -V, --version           Print program version.

Mandatory or optional arguments to long options are also mandatory or optional
for any corresponding short options.

GNU Astronomy Utilities home page: http://www.gnu.org/software/gnuastro/

Report bugs to bug-gnuastro@gnu.org.
EOF
}





# Output of `--version':
print_version() {
    cat <<EOF
$scriptname (GNU Astronomy Utilities) $version
Copyright (C) 2020-2022 Free Software Foundation, Inc.
License GPLv3+: GNU General public license version 3 or later.
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to the extent permitted by law.

Written/developed by Raul Infante-Sainz
EOF
}





# Output of `--cite':
print_citation() {
    empty="" # needed for the ascii art!
    cat <<EOF

Thank you for using $scriptname (GNU Astronomy Utilities) $version

Citations and acknowledgement are vital for the continued work on Gnuastro.

Please cite the following record(s) and add the acknowledgement statement below in your work to support us. Please note that different Gnuastro programs may have different corresponding papers. Hence, please check all the programs you used. Don't forget to also include the version as shown above for reproducibility.

Paper describing the creation of an extended PSF
------------------------------------------------
@ARTICLE{gnuastro-psf,
       author = {{Infante-Sainz}, Ra{\'u}l and {Trujillo}, Ignacio and {Rom{\'a}n}, Javier},
        title = "{The Sloan Digital Sky Survey extended point spread functions}",
      journal = {MNRAS},
         year = 2020,
        month = feb,
       volume = {491},
       number = {4},
        pages = {5317},
          doi = {10.1093/mnras/stz3111},
archivePrefix = {arXiv},
       eprint = {1911.01430},
 primaryClass = {astro-ph.IM},
       adsurl = {https://ui.adsabs.harvard.edu/abs/2020MNRAS.491.5317I},
      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
}

First paper introducing Gnuastro
--------------------------------
  @ARTICLE{gnuastro,
     author = {{Akhlaghi}, M. and {Ichikawa}, T.},
      title = "{Noise-based Detection and Segmentation of Nebulous Objects}",
    journal = {ApJS},
  archivePrefix = "arXiv",
     eprint = {1505.01664},
   primaryClass = "astro-ph.IM",
       year = 2015,
      month = sep,
     volume = 220,
        eid = {1},
      pages = {1},
        doi = {10.1088/0067-0049/220/1/1},
     adsurl = {https://ui.adsabs.harvard.edu/abs/2015ApJS..220....1A},
    adsnote = {Provided by the SAO/NASA Astrophysics Data System}
  }

Acknowledgement
---------------
This work was partly done using GNU Astronomy Utilities (Gnuastro,
ascl.net/1801.009) version $version. Work on Gnuastro has been funded by
the Japanese Ministry of Education, Culture, Sports, Science, and
Technology (MEXT) scholarship and its Grant-in-Aid for Scientific Research
(21244012, 24253003), the European Research Council (ERC) advanced grant
339659-MUSICOS, the Spanish Ministry of Economy and Competitiveness
(MINECO, grant number AYA2016-76219-P) and the NextGenerationEU grant
through the Recovery and Resilience Facility project
ICTS-MRR-2021-03-CEFCA.
                                               ,
                                              {|'--.
                                             {{\    \ $empty
      Many thanks from all                   |/\`'--./=.
      Gnuastro developers!                   \`\.---' \`\\
                                                  |\  ||
                                                  | |//
                                                   \//_/|
                                                   //\__/
                                                  //
                   (http://www.chris.com/ascii/) |/

EOF
}





# Functions to check option values and complain if necessary.
on_off_option_error() {
    if [ x"$2" = x ]; then
        echo "$scriptname: '$1' doesn't take any values"
    else
        echo "$scriptname: '$1' (or '$2') doesn't take any values"
    fi
    exit 1
}

check_v() {
    if [ x"$2" = x ]; then
        cat <<EOF
$scriptname: option '$1' requires an argument. Try '$scriptname --help' for more information
EOF
        exit 1;
    fi
}





# Separate command-line arguments from options. Then put the option
# value into the respective variable.
#
# OPTIONS WITH A VALUE:
#
#   Each option has three lines because we want to all common formats: for
#   long option names: `--longname value' and `--longname=value'. For short
#   option names we want `-l value', `-l=value' and `-lvalue' (where `-l'
#   is the short version of the hypothetical `--longname' option).
#
#   The first case (with a space between the name and value) is two
#   command-line arguments. So, we'll need to shift it two times. The
#   latter two cases are a single command-line argument, so we just need to
#   "shift" the counter by one. IMPORTANT NOTE: the ORDER OF THE LATTER TWO
#   cases matters: `-h*' should be checked only when we are sure that its
#   not `-h=*').
#
# OPTIONS WITH NO VALUE (ON-OFF OPTIONS)
#
#   For these, we just want the two forms of `--longname' or `-l'. Nothing
#   else. So if an equal sign is given we should definitely crash and also,
#   if a value is appended to the short format it should crash. So in the
#   second test for these (`-l*') will account for both the case where we
#   have an equal sign and where we don't.
inputs=""
while [ $# -gt 0 ]
do
    case "$1" in
        # Input parameters.
        -h|--hdu)            hdu="$2";                                  check_v "$1" "$hdu";  shift;shift;;
        -h=*|--hdu=*)        hdu="${1#*=}";                             check_v "$1" "$hdu";  shift;;
        -h*)                 hdu=$(echo "$1"  | sed -e's/-h//');        check_v "$1" "$hdu";  shift;;
        -n|--normradii)      normradii="$2";                            check_v "$1" "$normradii";  shift;shift;;
        -n=*|--normradii=*)  normradii="${1#*=}";                       check_v "$1" "$normradii";  shift;;
        -n*)                 normradii=$(echo "$1"  | sed -e's/-n//');  check_v "$1" "$normradii";  shift;;
        -W|--stampwidth)     stampwidth="$2";                           check_v "$1" "$stampwidth";  shift;shift;;
        -W=*|--stampwidth=*) stampwidth="${1#*=}";                      check_v "$1" "$stampwidth";  shift;;
        -W*)                 stampwidth=$(echo "$1"  | sed -e's/-W//'); check_v "$1" "$stampwidth";  shift;;
        -w|--corewidth)      corewidth="$2";                            check_v "$1" "$corewidth";  shift;shift;;
        -w=*|--corewidth=*)  corewidth="${1#*=}";                       check_v "$1" "$corewidth";  shift;;
        -w*)                 corewidth=$(echo "$1"  | sed -e's/-w//');  check_v "$1" "$corewidth";  shift;;
        -S|--segment)        segment="$2";                              check_v "$1" "$segment";  shift;shift;;
        -S=*|--segment=*)    segment="${1#*=}";                         check_v "$1" "$segment";  shift;;
        -S*)                 segment=$(echo "$1"  | sed -e's/-S//');    check_v "$1" "$segment";  shift;;
        -O|--mode)           mode="$2";                                 check_v "$1" "$mode";  shift;shift;;
        -O=*|--mode=*)       mode="${1#*=}";                            check_v "$1" "$mode";  shift;;
        -O*)                 mode=$(echo "$1"  | sed -e's/-O//');       check_v "$1" "$mode";  shift;;
        -c|--center)         center="$2";                               check_v "$1" "$center";  shift;shift;;
        -c=*|--center=*)     center="${1#*=}";                          check_v "$1" "$center";  shift;;
        -c*)                 center=$(echo "$1"  | sed -e's/-c//');     check_v "$1" "$center";  shift;;
        -N|--normop)         normop="$2";                               check_v "$1" "$normop";  shift;shift;;
        -N=*|--normop=*)     normop="${1#*=}";                          check_v "$1" "$normop";  shift;;
        -N*)                 normop=$(echo "$1"  | sed -e's/-N//');     check_v "$1" "$normop";  shift;;
        -s|--sigmaclip)      sigmaclip="$2";                            check_v "$1" "$sigmaclip";  shift;shift;;
        -s=*|--sigmaclip=*)  sigmaclip="${1#*=}";                       check_v "$1" "$sigmaclip";  shift;;
        -s*)                 sigmaclip=$(echo "$1"  | sed -e's/-s//');  check_v "$1" "$sigmaclip";  shift;;
        -Q|--axisratio)      axisratio="$2";                            check_v "$1" "$axisratio";  shift;shift;;
        -Q=*|--axisratio=*)  axisratio="${1#*=}";                       check_v "$1" "$axisratio";  shift;;
        -Q*)                 axisratio=$(echo "$1"  | sed -e's/-Q//');  check_v "$1" "$axisratio";  shift;;
        -p|--positionangle)     positionangle="$2";                            check_v "$1" "$positionangle";  shift;shift;;
        -p=*|--positionangle=*) positionangle="${1#*=}";                       check_v "$1" "$positionangle";  shift;;
        -p*)                    positionangle=$(echo "$1"  | sed -e's/-p//');  check_v "$1" "$positionangle";  shift;;


        # Output parameters
        -k|--keeptmp)     keeptmp=1; shift;;
        -k*|--keeptmp=*)  on_off_option_error --keeptmp -k;;
        -t|--tmpdir)      tmpdir="$2";                          check_v "$1" "$tmpdir";  shift;shift;;
        -t=*|--tmpdir=*)  tmpdir="${1#*=}";                     check_v "$1" "$tmpdir";  shift;;
        -t*)              tmpdir=$(echo "$1" | sed -e's/-t//'); check_v "$1" "$tmpdir";  shift;;
        -o|--output)      output="$2";                          check_v "$1" "$output"; shift;shift;;
        -o=*|--output=*)  output="${1#*=}";                     check_v "$1" "$output"; shift;;
        -o*)              output=$(echo "$1" | sed -e's/-o//'); check_v "$1" "$output"; shift;;

        # Non-operating options.
        -q|--quiet)       quiet="--quiet"; shift;;
        -q*|--quiet=*)    on_off_option_error --quiet -q;;
        -?|--help)        print_help; exit 0;;
        -'?'*|--help=*)   on_off_option_error --help -?;;
        -V|--version)     print_version; exit 0;;
        -V*|--version=*)  on_off_option_error --version -V;;
        --cite)           print_citation; exit 0;;
        --cite=*)         on_off_option_error --cite;;

        # Unrecognized option:
        -*) echo "$scriptname: unknown option '$1'"; exit 1;;

        # Not an option (not starting with a `-'): assumed to be input FITS
        # file name.
        *) if [ x"$inputs" = x ]; then inputs="$1"; else inputs="$inputs $1"; fi; shift;;
    esac

done





# Basic sanity checks
# ===================

# If an input image is not given at all.
if [ x"$inputs" = x ]; then
    cat <<EOF
$scriptname: no input FITS image files. Run with '--help' for more information on how to run
EOF
    exit 1
elif [ ! -f "$inputs" ]; then
    echo "$scriptname: $inputs: No such file or directory"
    exit 1
fi

# If a stamp width (--stampwidth) is not given at all.
if [ x"$stampwidth" = x ]; then
    cat <<EOF
$scriptname: no stamp width provided. A stamp size (width) has to be specified with '--stampwidth' (or '-w')
EOF
    exit 1
fi

# If center coordinates (--center) is not given at all.
if [ x"$center" = x ]; then
    cat <<EOF
$scriptname: no center coordinates provided. Please use '--center' ('-c') to give the central coordinates
EOF
    exit 1
else
    ncenter=$(echo $center | awk 'BEGIN{FS=","}END{print NF}')
    if [ x$ncenter != x2 ]; then
        cat <<EOF
$scriptname: '--center' (or '-c') only take two values, but $ncenter were given in '$center'
EOF
        exit 1
    fi
fi

# If a normalization range is not given at all.
if [ x"$normradii" = x ]; then
    cat <<EOF
$scriptname: warning: no ring of normalization provided. the stamp won't be normalized
EOF
else
    nnormradii=$(echo $normradii | awk 'BEGIN{FS=","}END{print NF}')
    if [ x$nnormradii != x2 ]; then
        cat <<EOF
$scriptname: '--normradii' (or '-n') only take two values, but $nnormradii were given in '$normradii'
EOF
        exit 1
    fi
fi

# If mode (--mode) is not given at all.
if [ x"$mode" = x ]; then
    cat <<EOF
$scriptname: no coordinate mode provided. The '--mode' ('-O') takes one of the following two values: 'img' (for pixel coordinates) or 'wcs' (for celestial coordinates)
EOF
    exit 1

# Make sure the value to '--mode' is either 'wcs' or 'img'. Note: '-o'
# means "or" and is preferred to '[ ] || [ ]' because only a single
# invocation of 'test' is done. Run 'man test' for more.
elif [ "$mode" = wcs     -o      $mode = "img" ]; then
    junk=1
else
    cat <<EOF
$scriptname: value to '--mode' (or '-O') is not recognized ('$mode'). This option takes one of the following two values: 'img' (for pixel coordinates) or 'wcs' (for celestial coordinates)
EOF
    exit 1
fi





# Basic parameters: coordinates and normalization radii
# -----------------------------------------------------
#
# Obtain the coordinates of the center as well as the normalization radii
# ring from the command line arguments.
xcoord=$(echo "$center" | awk 'BEGIN{FS=","} {print $1}')
ycoord=$(echo "$center" | awk 'BEGIN{FS=","} {print $2}')

normradiusmin=$(echo "$normradii" | awk 'BEGIN{FS=","} {print $1}')
normradiusmax=$(echo "$normradii" | awk 'BEGIN{FS=","} {print $2}')

# With the center coordinates, generate a specific label for the object
# consisting in its coordinates.
objectid="$xcoord"_"$ycoord"





# Define a temporal directory and thefinal output file
# ----------------------------------------------------
#
# Construct the temporary directory. If the user does not specify any
# directory, then a default one with the base name of the input image will
# be constructed.  If the user set the directory, then make it. This
# directory will be deleted at the end of the script if the user does not
# want to keep it (with the `--keeptmp' option).

# The final output stamp is also defined here if the user does not provide
# an explicit name. If the user has defined a specific path/name for the
# output, it will be used for saving the output file. If the user does not
# specify a output name, then a default value containing the center and
# mode will be generated.
bname_prefix=$(basename "$inputs" | sed 's/\.fits/ /' | awk '{print $1}')
if [ x"$tmpdir" = x ]; then \
  tmpdir=$(pwd)/"$bname_prefix"_stamp
fi

if [ -d "$tmpdir" ]; then
  junk=1
else
  mkdir -p "$tmpdir"
fi

# Output
if [ x"$output" = x ]; then
  output="$bname_prefix"_stamp_$objectid.fits
fi





# Transform WCS to IMG center coordinates
# ---------------------------------------
#
# If the original coordinates have been given in WCS or celestial units
# (RA/DEC), then transform them to IMG (pixel). Here, this is done by using
# the WCS information from the original input image. If the original
# coordinates were done in IMG, then just use them.
if [ "$mode" = wcs ]; then
  xycenter=$(echo "$xcoord,$ycoord" \
                  | asttable  --column='arith $1 $2 wcs-to-img' \
                              --wcsfile=$inputs --wcshdu=$hdu $quiet)
  xcenter=$(echo "$xycenter" | awk '{print $1}')
  ycenter=$(echo "$xycenter" | awk '{print $2}')
else
  xcenter=$xcoord
  ycenter=$ycoord
fi





# Crop the original image around the object
# -----------------------------------------
#
# Crop the object around its center with the given stamp size width.
cropped=$tmpdir/cropped-$objectid.fits
astcrop $inputs --hdu=$hdu --mode=img \
                --center=$xcenter,$ycenter \
                --width=$stampwidth --output=$cropped $quiet





# Function to find label of desired region
# ----------------------------------------
#
# Given a certain extension ('CLUMPS' or 'OBJECTS'), find the respective
# label.
find_central_label() {

  # Input arguments
  hdu=$1

  # Set the parameters.
  case $hdu in
      CLUMPS) labname=clump;;
      OBJECTS) labname=object;;
      *) cat <<EOF
$scriptname: ERROR: a bug! Please contact us at bug-gnuastro@gnu.org to fix the problem. The first argument to 'find_central_label' is not 'CLUMPS' or 'OBJECTS'
EOF
         exit 1;;
  esac

  # Crop the labeled image.
  cropped_core=$tmpdir/cropped-core-$objectid.fits
  astcrop $segment --hdu=$hdu --mode=img $quiet --width=$corewidth \
          --center=$xcenter,$ycenter --output=$cropped_core
  lab=$(astarithmetic $cropped_core unique --quiet)
  numlab=$(echo "$lab" | awk '{print NF}')
  if [ $numlab != 1 ]; then
    cat <<EOF
$scriptname: ERROR: there is more than one $labname label in the core region around the given coordinate (specified by '--corewidth=$corewidth', containing $labname labels: $lab). Therefore it is not possible to unambiguously identify a single $labname. Please decrease the box size given to '--corewidth'.
EOF
    exit 1
  fi

  # Write the output in a file.
  echo $lab > $tmpdir/cropped-core-$labname-$objectid.txt
}




# Crop and unlabel the segmentation image
# ---------------------------------------
#
# If the user provides a segmentation image, treat it appropiately in order
# to mask all objects that are not the central one. If not, just consider
# that the cropped and masked image is the cropped (not masked) image. The
# process is as follow:
#   - Crop the original mask image.
#   - Crop the original mask image to a central region (core), in order to
#     compute what is the central object id. This is necessary to unmask
#     this object.
#   - Compute what is the central object value, using the median value.
#   - In the original cropped mask, convert all pixels belonging to the
#     central object to zeros. By doing this, the central object becomes as
#     sky.
#   - Mask all non zero pixels in the mask image as nan values.
if [ x"$segment" != x ]; then

  # Find the object and clump labels of the target.
  find_central_label CLUMPS
  find_central_label OBJECTS
  clab=$(cat $tmpdir/cropped-core-clump-$objectid.txt)
  olab=$(cat $tmpdir/cropped-core-object-$objectid.txt)

  # Crop the object and clump image to same size as desired stamp.
  cropclp=$tmpdir/cropped-clumps-$objectid.fits
  cropobj=$tmpdir/cropped-objects-$objectid.fits
  astcrop $segment --hdu=OBJECTS --mode=img \
          --center=$xcenter,$ycenter \
          --width=$stampwidth --output=$cropobj $quiet
  astcrop $segment --hdu=CLUMPS --mode=img \
          --center=$xcenter,$ycenter \
          --width=$stampwidth --output=$cropclp $quiet

  # Mask all the undesired regions.
  cropped_masked=$tmpdir/cropped-masked-$objectid.fits
  astarithmetic $cropped --hdu=1 set-i --output=$cropped_masked \
                $cropobj --hdu=1 set-o \
                $cropclp --hdu=1 set-c \
                                       \
                c o $olab ne 0 where c $clab eq -1 where 0 gt set-cmask \
                o o $olab eq 0 where set-omask \
                i omask cmask or nan where
else
  cropped_masked=$cropped
fi





# Compute the radial profile and the normalization value
# ------------------------------------------------------
#
# Only if the the user has specified a ring of normalization (--normradii).
# Otherwise set the normalization value equal to 1.0 (no normalization).
if [ x"$normradiusmin" != x   -a   x"$normradiusmax" != x ]; then

    # Generate the radial profile of the stamp, since it has been already
    # centered on the center of the object, it is not necessary to give the
    # center coordinates. If the user specifies a maximum radius, use it.
    # Otherwise, compute the radial profile up to the outer part of the
    # ring for the normalization (to not wast CPU time). If the user
    # specifies sigma clip parameters, use them.
    radialprofile=$tmpdir/rprofile-$objectid.fits
    maxr=$(echo "$normradiusmax" | awk '{print $1+1}')
    if [ x"$sigmaclip" = x ]; then finalsigmaclip=""
    else                           finalsigmaclip="--sigmaclip=$sigmaclip";
    fi
    astscript-radial-profile $cropped_masked --hdu=1 --rmax=$maxr \
                             --measure=$normop $finalsigmaclip \
                             --positionangle=$positionangle \
                             --tmpdir=$tmpdir --keeptmp \
                             --axisratio=$axisratio \
                             --output=$radialprofile $quiet

    # The normalization value is computed from the radial profile in between
    # the two radius specified by the user. In this case, the option to give
    # sigmaclip parameters to 'aststatistics' is different.
    if [ x"$sigmaclip" = x ]; then finalsigmaclip=""
    else                           finalsigmaclip="--sclipparams=$sigmaclip";
    fi
    normvalue=$(asttable $radialprofile $quiet \
                         --range=1,$normradiusmin,$normradiusmax \
                         | aststatistics --column=2 --$normop \
                                         $finalsigmaclip -q)
else
    normvalue=1.0
fi





# Normalize the stamp
# -------------------
#
# Before applying the normalization factor to the masked stamps, we should
# make sure the the normalization value has a 'float32' type. This is
# because Arithmetic will interpret its type based on the number of digits
# it has. Therefore something like 1.7204 will be read as float32 (1.7204),
# but 1.868915 will be read as 'float64'. In the latter case, the output
# will become 'float64' also, which will cause problems later when we want
# to stack them together (some images will be 'float32', some 'float64').
astarithmetic $cropped_masked --hdu=1 $normvalue float32 / \
              --output=$output $quiet





# Remove temporary files
# ----------------------
#
# If the user does not specify to keep the temporal files with the option
# `--keeptmp', then remove the whole directory.
if [ $keeptmp = 0 ]; then
    rm -r $tmpdir
fi
