import sys
import os
import time
from signal import signal, SIGINT
import requests
import socket
import warnings

import numpy as np

import matplotlib.pylab as plt
from matplotlib.patches import Circle

from astropy.time import Time
from astropy.coordinates import Angle
from astropy import units as u

from trm import cline
from trm.cline import Cline

import hipercam as hcam
from hipercam import defect, fringe, mpl, spooler
from hipercam.mpl import Params

try:
    import lacosmic
except ModuleNotFoundError:
    lacosmic = None

# colour for setup windows. works for me at least
# but may need input from Stu
COL_SETUP = (0.8, 0.0, 0.0)

__all__ = [
    "nrtplot",
]

######################################
#
# nrtplot -- display of multiple images
#
######################################


def nrtplot(args=None):
    """``nrtplot [source] (run first [twait tmax] | flist) trim ([ncol
    nrow]) (ccd (nx)) [imwidth pause tupdate plotall] bias [lowlevel
    highlevel] dark flat fmap (fpair [nhalf rmin rmax]) defect setup
    [drurl cmap imwidth imheight memory] msub iset (ilo ihi | plo phi)
    xlo xhi ylo yhi profit [method beta fwhm fwhm_min shbox smooth
    fhbox hmin read gain thresh (fwnmax fwymax fwwidth fwheight)]
    [fcosmic] (contrast crthresh nbthresh sathresh readout)``

    This is 'nrtplot' "new" rtplot, a matplotlib-based replacement for
    the current PGPLOT one. Under development.

    Plots a sequence of images as a movie in near 'real time', hence
    'rt'. Designed to be used to look at images coming in while at the
    telescope, 'nrtplot' comes with many options, a large number of
    which are hidden by default, and many of which are only prompted
    if other arguments are set correctly. If you want to see them all,
    invoke as 'nrtplot prompt'.  This is worth doing once to know
    nrtplot's capabilities.

    nrtplot can source data from both the ULTRACAM and HiPERCAM
    servers, from local 'raw' ULTRACAM and HiPERCAM files (i.e. .xml +
    .dat for ULTRACAM, 3D FITS files for HiPERCAM) and from lists of
    HiPERCAM '.hcm' files.

    nrtplot optionally allows the selection of targets to be fitted
    with gaussian or moffat profiles, and, if successful, will plot
    circles of 2x the measured FWHM in green over the selected
    targets. In this case it can also plot the history of the FWHMs as
    an aid to focussing. You can also click on targets on the fly to
    get one-off measurement which are just printed to the screen. See
    below for a note on this.

    You can re-size and zoom & pan during plotting, although it takes
    some getting used owing to the way the program blocks. For re-sizing,
    you may need to keep the left-button of the mouse down for a while
    for anything to happen. See the parameter 'memory' below which can
    help.

    When you first run it, try to get the initial fit parameters about
    right; it doesn't like to be miles out at the start for some reason
    and will refuse to fit.

    Experimental: I have added options to apply L.A.cosmic just to see what
    it returns at the moment.

    Parameters:

        source : str [hidden]
           Data source, five options:

             |  'hs' : HiPERCAM server
             |  'hl' : local HiPERCAM FITS file
             |  'us' : ULTRACAM server
             |  'ul' : local ULTRACAM .xml/.dat files
             |  'hf' : list of HiPERCAM hcm FITS-format files

           'hf' is used to look at sets of frames generated by 'grab'
           or converted from foreign data formats. The standard
           start-off default for ``source`` can be set using the
           environment variable HIPERCAM_DEFAULT_SOURCE. e.g. in bash
           :code:`export HIPERCAM_DEFAULT_SOURCE="us"` would ensure it
           always started with the ULTRACAM server by default. If
           unspecified, it defaults to 'hl'.

        run : str [if source ends 's' or 'l']
           run number to access, e.g. 'run034'

        flist : str [if source ends 'f']
           name of file list

        first : int [if source ends 's' or 'l']
           exposure number to start from. 1 = first frame; set = 0 to always
           try to get the most recent frame (if it has changed).  For data
           from the |hiper| server, a negative number tries to get a frame not
           quite at the end.  i.e. -10 will try to get 10 from the last
           frame. This is mainly to sidestep a difficult bug with the
           acquisition system.

        twait : float [if source ends 's' or 'l'; hidden]
           time to wait between attempts to find a new exposure, seconds.

        tmax : float [if source ends 's' or 'l'; hidden]
           maximum time to wait between attempts to find a new exposure,
           seconds.

        trim : bool [if source starts with 'u']
           True to trim columns and/or rows off the edges of windows nearest
           the readout which can sometimes contain bad data.

        ncol : int [if trim, hidden]
           Number of columns to remove (on left of left-hand window, and right
           of right-hand windows)

        nrow : int [if trim, hidden]
           Number of rows to remove (bottom of windows)

        ccd : str
           CCD(s) to plot, '0' for all, '1 3' to plot '1' and '3' only, etc.

        nx : int [if more than 1 CCD]
           number of panels across to display for the image display.

        pause : float [hidden]
           seconds to pause between frames (defaults to 0). This is in addition
           to any time taken per frame to read and display it, so is just to slow
           things not to set an exact rate.

        tupdate : float [hidden]
           target number of seconds between plot updates. This is to
           make the plots more responsive and get round issues with
           plots freezing until the next frame comes through which can
           be especially annoying with virtual desktops. Experiments
           with ULTRACAM@NTT indicated 0.2 to 0.6 secs or so works OK.
           Main impact seems to be on drpc CPU.

        plotall : bool [hidden]
           plot all frames regardless of status (i.e. including blank frames
           when nskips are enabled (defaults to False). The profile fitting
           will still be disabled for bad frames.

        bias : str
           Name of bias frame to subtract, 'none' to ignore.

        lowlevel : float [hidden]
           Level below which a warning about low bias levels is warned. Set=0
           to ignore. Applied to first window of first CCD. 2000 about
           right for ULTRACAM.

        highlevel : float [hidden]
           Level above which a warning about high bias levels is warned. Set=0
           to ignore. Applied to first window of first CCD. 3500 about
           right for ULTRACAM.

        dark : str
           Name of dark frame to subtract, 'none' to ignore

        flat : str
           Name of flat field to divide by, 'none' to ignore. Should normally
           only be used in conjunction with a bias, although it does allow you
           to specify a flat even if you haven't specified a bias.

        fmap : str
           Name of fringe map (see e.g. `makefringe`), 'none' to ignore.

        fpair : str [if fmap is not 'none']
           Name of fringe pair file (see e.g. `setfringe`). Required if
           a fringe map has been specified.

        nhalf : int [if fmap is not 'none', hidden]
           When calculating the differences for fringe measurement,
           a region extending +/-nhalf binned pixels will be used when
           measuring the amplitudes. Basically helps the stats.

        rmin : float [if fmap is not 'none', hidden]
           Minimum individual ratio to accept prior to calculating the overall
           median in order to reduce the effect of outliers. Although all ratios
           should be positive, you might want to set this a little below zero
           to allow for some statistical fluctuation.

        rmax : float [if fmap is not 'none', hidden]
           Maximum individual ratio to accept prior to calculating the overall
           median in order to reduce the effect of outliers. Probably typically
           < 1 if fringe map was created from longer exposure data.

        defect : str
           Name of defect file, 'none' to ignore.

        setup : bool
           True/yes to access the current windows from hdriver. Useful
           during observing when setting up windows, but not normally
           otherwise. Next argument (hidden) is the URL to get to
           hdriver. Once setup, you should probably turn this off to
           avoid overloading hdriver, especially if in drift mode as
           it makes a request for the windows for every frame.

        drurl : str [if setup; hidden]
           URL needed to access window setting from the camera
           driver (ultracam, ultraspec, hipercam). The internal server
           in the camera driver must be switched on which can be done
           from the GUI.

        cmap : str [hidden]
           The matplotlib colour map to use. "Greys" gives the usual greyscale.
           "none" will give whatever the current default is. Many other choices:
           "viridis", "jet", "hot", "Oranges", etc. Enter an invalid one and
           the program will fail but return a huge list of possibles in the
           process. Note that some will not work well with the colours used
           to plot annotations.

        imwidth : float [hidden]
           image display plot width in inches (0 for default, and note
           that this will also make the height go to its default)

        imheight : float [hidden]
           image display plot height in inches (0 for default, and note
           that this will also make the width go to its default)

        memory : bool [hidden]
           stores the image display dimensions at the end as defaults for
           imwidth and imheight next time (or not). This helps when you are
           repeatedly re-starting nrtplot for the same sort of display.

        msub : bool
           subtract the median from each window before scaling for the
           image display or not. This happens after any bias subtraction.

        iset : str [single character]
           determines how the intensities are determined. There are three
           options: 'a' for automatic simply scales from the minimum to the
           maximum value found on a per CCD basis. 'd' for direct just takes
           two numbers from the user. 'p' for percentile dtermines levels
           based upon percentiles determined from the entire CCD on a per CCD
           basis.

        ilo : float [if iset='d']
           lower intensity level

        ihi : float [if iset='d']
           upper intensity level

        plo : float [if iset='p']
           lower percentile level

        phi : float [if iset='p']
           upper percentile level

        xlo : float
           left-hand X-limit for plot, initially at least since it is
           possible to re-size. For iset='p' these limits also set the
           region of the frame over which the percentil will be
           calculated. You will usually want yhi-ylo ~ xhi-xlo in
           magnitude because the aspect ratio is preserved.

        xhi : float
           right-hand X-limit for plot (can be < xlo to invert the display)

        ylo : float
           lower Y-limit for plot

        yhi : float
           upper Y-limit for plot (can be < ylo)

        profit : bool
           If profit=True, you say yes, then you will get the option
           of selecting a fixed set of objects with a cursor, and the
           program will then attempt to track these from frame to
           frame, and fit their profile. You may need to adjust
           'first' to see anything. Depending on 'fnmax', it will also
           then plot a history of the FWHM measurements, which is
           useful for focussing. The parameters used for profile fits
           are hidden and you may want to invoke the command with
           'prompt' the first time you try profile fitting. NB There
           is also an option for "on-the-fly" fits of any target: see
           the note below.

        method : str [hidden]
           this defines the profile fitting method, either a gaussian or a
           moffat profile. The latter is usually best.

        beta : float [method == 'm'; hidden]
           default Moffat exponent

        fwhm : float [hidden]
           default FWHM, unbinned pixels. Do try to get this about right as
           it affects whether the profile can be fitted at all.

        fwhm_min : float [hidden]
           minimum FWHM to allow, unbinned pixels.

        shbox : float [hidden]
           half width of box for searching for a star, unbinned
           pixels. The above-threshold target closest to the centre of
           the box in a region +/- shbox around an intial position
           will be selected. It may not be the brightest, depending
           upon your threshold settings, so use those to filter faint
           objects. If profit=True, 'shbox' should be large enough to
           allow for likely changes in position from frame to frame,
           but not too large to avoid jumping to brighter targets or
           possibly cosmic rays.

        smooth : float [hidden]
           FWHM for gaussian smoothing, binned pixels. The initial position
           for fitting is determined by finding the maximum flux in a smoothed
           version of the image in a box of width +/- shbox around the starter
           position. Typically should be comparable to the stellar width. Its
           main purpose is to combat cosmic rays which tend only to occupy a
           single pixel.

        fhbox : float [hidden]
           half width of box for profile fit, unbinned pixels. The fit box is
           centred on the position located by the initial search. It should
           normally be > ~2x the expected FWHM, and usually smaller than shbox

        hmin : float [hidden]
           height threshold to accept a fit. If the height is below this
           value, the position will not be updated. This is to help in cloudy
           conditions. The limit is applied to the image after it has been
           smoothed to make less vulnerable to seeing fluctuations. This
           can mean it can be quite small.

        read : float [hidden]
           readout noise, RMS ADU, for assigning uncertainties. Set to a -ve
           value to try to ascertain on the fly; this is advisable either if
           you don't have a bias or you apply 'msub', and it will default to
           this if not specified in this case. The value returned in this case
           includes sky noise, i.e it should be roughly sqrt(R**2+S/G) where
           R is the true read noise, S are the sky counts per pixel, and G the
           gain. If read is set -ve, two fits are carried out per target. The
           second of these should usually be pretty fast. The first is carried
           out with an assumed large read noise of 20 in order to soften the
           weights. The results reported apply to the second fit.

        gain : float [hidden]
           gain, ADU/count, for assigning uncertainties.

        thresh : float [hidden]
           sigma rejection threshold for fits

        fwnmax : int [if profit; hidden]
           maximum number of frames to buffer the FWHM plot (< 2 ==>
           no plot). It allows the FWHMs of multiple stars to be
           tracked which might be useful for focussing. The plot
           starts to "slide" once the buffers are filled.

        fwymax : float [if profit; hidden]
           maximum FWHM to plot on Y-axis [unbinned pixels]. It is
           possible to alter the value of the fly by zooming and panning
           the plot.

        fwwidth : float [if profit; hidden]
           FWHM display plot width in inches (0 for default, which will
           also cause the height to go to its default value)

        fwheight : float [if profit; hidden]
           FWHM display plot height in inches (0 for default, which will
           also cause the width to go to its default value)

        fcosmic : bool [hidden, defaults to False]
           True to identify cosmic rays using L.A.cosmic (van Dokkum 2001,
           PASP, 113, 1420) as implemented in Python package "lascosmic".
           Experimental at the moment. Just returns the number of cosmic
           rays found.

        crthresh : float
           Threshold signal-to-noise ratio of cosmic rays

        nbthresh : float
           Threshold signal-to-noise ratio of pixels neighbouring
           cosmic rays

        sathresh : float
           Threshold to flag saturated pixels

        readout : float
           readout noise, RMS ADU, for assigning uncertainties. (NB
           at the moment this is separate from "read", but they may
           be merged at some point)

    .. Note::

        To help with *dithered* long exposures especially, for which
        'profit' performs very poorly, clicking on any object will
        attempt a one-off profile fit, with results that are reported
        to the terminal. It's a little cludgy in that the profile fits
        are only carried out and reported *after* the next frame has
        been plotted, but it's better than nothing.

        If you don't apply bias subtraction and/or you do subtract the
        median, then the usual method to estimate the read noise fails.
        In that case, setting 'read' < 0 will force an attempt to measure
        the read noise from fluctuations in the sky. You need to ensure
        that fhbox is large enough to encompass sky regions well away from
        the target in this case.

    """

    command, args = cline.script_args(args)

    # get the inputs
    with Cline("HIPERCAM_ENV", ".hipercam", command, args) as cl:

        # register parameters
        cl.register("source", Cline.GLOBAL, Cline.HIDE)
        cl.register("run", Cline.GLOBAL, Cline.PROMPT)
        cl.register("first", Cline.LOCAL, Cline.PROMPT)
        cl.register("trim", Cline.GLOBAL, Cline.PROMPT)
        cl.register("ncol", Cline.GLOBAL, Cline.HIDE)
        cl.register("nrow", Cline.GLOBAL, Cline.HIDE)
        cl.register("twait", Cline.LOCAL, Cline.HIDE)
        cl.register("tmax", Cline.LOCAL, Cline.HIDE)
        cl.register("flist", Cline.LOCAL, Cline.PROMPT)
        cl.register("ccd", Cline.LOCAL, Cline.PROMPT)
        cl.register("nx", Cline.LOCAL, Cline.PROMPT)
        cl.register("pause", Cline.LOCAL, Cline.HIDE)
        cl.register("tupdate", Cline.LOCAL, Cline.HIDE)
        cl.register("plotall", Cline.LOCAL, Cline.HIDE)
        cl.register("bias", Cline.GLOBAL, Cline.PROMPT)
        cl.register("lowlevel", Cline.GLOBAL, Cline.HIDE)
        cl.register("highlevel", Cline.GLOBAL, Cline.HIDE)
        cl.register("dark", Cline.GLOBAL, Cline.PROMPT)
        cl.register("flat", Cline.GLOBAL, Cline.PROMPT)
        cl.register("fmap", Cline.GLOBAL, Cline.PROMPT)
        cl.register("fpair", Cline.GLOBAL, Cline.PROMPT)
        cl.register("nhalf", Cline.GLOBAL, Cline.HIDE)
        cl.register("rmin", Cline.GLOBAL, Cline.HIDE)
        cl.register("rmax", Cline.GLOBAL, Cline.HIDE)
        cl.register("defect", Cline.GLOBAL, Cline.PROMPT)
        cl.register("setup", Cline.GLOBAL, Cline.PROMPT)
        cl.register("drurl", Cline.GLOBAL, Cline.HIDE)
        cl.register("cmap", Cline.LOCAL, Cline.HIDE)
        cl.register("imwidth", Cline.LOCAL, Cline.HIDE)
        cl.register("imheight", Cline.LOCAL, Cline.HIDE)
        cl.register("memory", Cline.LOCAL, Cline.HIDE)
        cl.register("msub", Cline.GLOBAL, Cline.PROMPT)
        cl.register("iset", Cline.GLOBAL, Cline.PROMPT)
        cl.register("ilo", Cline.GLOBAL, Cline.PROMPT)
        cl.register("ihi", Cline.GLOBAL, Cline.PROMPT)
        cl.register("plo", Cline.GLOBAL, Cline.PROMPT)
        cl.register("phi", Cline.LOCAL, Cline.PROMPT)
        cl.register("xlo", Cline.GLOBAL, Cline.PROMPT)
        cl.register("xhi", Cline.GLOBAL, Cline.PROMPT)
        cl.register("ylo", Cline.GLOBAL, Cline.PROMPT)
        cl.register("yhi", Cline.GLOBAL, Cline.PROMPT)
        cl.register("profit", Cline.LOCAL, Cline.PROMPT)
        #        cl.register("nxf", Cline.LOCAL, Cline.HIDE)
        cl.register("method", Cline.LOCAL, Cline.HIDE)
        cl.register("beta", Cline.LOCAL, Cline.HIDE)
        cl.register("fwhm", Cline.LOCAL, Cline.HIDE)
        cl.register("fwhm_min", Cline.LOCAL, Cline.HIDE)
        cl.register("shbox", Cline.LOCAL, Cline.HIDE)
        cl.register("smooth", Cline.LOCAL, Cline.HIDE)
        cl.register("fhbox", Cline.LOCAL, Cline.HIDE)
        cl.register("hmin", Cline.LOCAL, Cline.HIDE)
        cl.register("read", Cline.LOCAL, Cline.HIDE)
        cl.register("gain", Cline.LOCAL, Cline.HIDE)
        cl.register("thresh", Cline.LOCAL, Cline.HIDE)
        cl.register("fwnmax", Cline.LOCAL, Cline.HIDE)
        cl.register("fwymax", Cline.LOCAL, Cline.HIDE)
        cl.register("fwwidth", Cline.LOCAL, Cline.HIDE)
        cl.register("fwheight", Cline.LOCAL, Cline.HIDE)
        cl.register("fcosmic", Cline.LOCAL, Cline.HIDE)
        cl.register("crthresh", Cline.LOCAL, Cline.PROMPT)
        cl.register("nrthresh", Cline.LOCAL, Cline.PROMPT)
        cl.register("sathresh", Cline.LOCAL, Cline.PROMPT)
        cl.register("readout", Cline.LOCAL, Cline.PROMPT)

        # get inputs
        default_source = os.environ.get("HIPERCAM_DEFAULT_SOURCE", "hl")
        source = cl.get_value(
            "source",
            "data source [hs, hl, us, ul, hf]",
            default_source,
            lvals=("hs", "hl", "us", "ul", "hf"),
        )

        # set some flags
        server_or_local = source.endswith("s") or source.endswith("l")

        if server_or_local:
            resource = cl.get_value("run", "run name", "run005")
            if source == "hs":
                first = cl.get_value("first", "first frame to plot", 1)
            else:
                first = cl.get_value("first", "first frame to plot", 1, 0)

            twait = cl.get_value(
                "twait", "time to wait for a new frame [secs]", 1.0, 0.0
            )
            tmax = cl.get_value(
                "tmax", "maximum time to wait for a new frame [secs]", 10.0, 0.0
            )

        else:
            resource = cl.get_value(
                "flist", "file list", cline.Fname("files.lis", hcam.LIST)
            )
            first = 1

        trim = cl.get_value("trim", "do you want to trim edges of windows?", True)
        if trim:
            ncol = cl.get_value("ncol", "number of columns to trim from windows", 0)
            nrow = cl.get_value("nrow", "number of rows to trim from windows", 0)
        else:
            ncol, nrow = None, None

        # define the panel grid. first get the labels and maximum dimensions
        ccdinf = spooler.get_ccd_pars(source, resource)

        nxdef = cl.get_default("nx", 3)

        if len(ccdinf) > 1:
            ccd = cl.get_value("ccd", "CCD(s) to plot [0 for all]", "0")
            if ccd == "0":
                ccds = list(ccdinf.keys())
            else:
                ccds = ccd.split()
                check = set(ccdinf.keys())
                if not set(ccds) <= check:
                    raise hcam.HipercamError("At least one invalid CCD label supplied")

            if len(ccds) > 1:
                nxdef = min(len(ccds), nxdef)
                cl.set_default("nx", nxdef)
                nx = cl.get_value("nx", "number of panels in X", 3, 1)
            else:
                nx = 1
        else:
            nx = 1
            ccds = list(ccdinf.keys())

        cl.set_default("pause", 0.0)
        pause = cl.get_value(
            "pause", "time delay to add between" " frame plots [secs]", 0.0, 0.0
        )

        tupdate = cl.get_value(
            "tupdate", "target time between plot updates [secs]", 0.2, 0.01
        )

        cl.set_default("plotall", False)
        plotall = cl.get_value(
            "plotall", "plot all frames," " regardless of status?", False
        )

        # bias frame (if any)
        bias = cl.get_value(
            "bias",
            "bias frame ['none' to ignore]",
            cline.Fname("bias", hcam.HCAM),
            ignore="none",
        )
        if bias is not None:
            # read the bias frame
            bias = hcam.MCCD.read(bias)
            fprompt = "flat frame ['none' to ignore]"
        else:
            fprompt = "flat frame ['none' is normal choice with no bias]"

        lowlevel = cl.get_value(
            "lowlevel", "bias level lower limit for warnings", 2000.0
        )

        highlevel = cl.get_value(
            "highlevel", "bias level upper limit for warnings", 3500.0
        )

        # dark (if any)
        dark = cl.get_value(
            "dark",
            "dark frame to subtract ['none' to ignore]",
            cline.Fname("dark", hcam.HCAM),
            ignore="none",
        )
        if dark is not None:
            # read the dark frame
            dark = hcam.MCCD.read(dark)

        # flat (if any)
        flat = cl.get_value(
            "flat", fprompt, cline.Fname("flat", hcam.HCAM), ignore="none"
        )
        if flat is not None:
            # read the flat frame
            flat = hcam.MCCD.read(flat)

        # fringe file (if any)
        fmap = cl.get_value(
            "fmap",
            "fringe map ['none' to ignore]",
            cline.Fname("fmap", hcam.HCAM),
            ignore="none",
        )
        if fmap is not None:
            # read the fringe map
            fmap = hcam.MCCD.read(fmap)
            fpair = cl.get_value(
                "fpair", "fringe pair file", cline.Fname("fpair", hcam.FRNG)
            )
            fpair = fringe.MccdFringePair.read(fpair)

            nhalf = cl.get_value(
                "nhalf", "half-size of fringe measurement regions", 2, 0
            )
            rmin = cl.get_value("rmin", "minimum fringe pair ratio", -0.2)
            rmax = cl.get_value("rmax", "maximum fringe pair ratio", 1.0)

        # defect file (if any)
        dfct = cl.get_value(
            "defect",
            "defect file ['none' to ignore]",
            cline.Fname("defect", hcam.DFCT),
            ignore="none",
        )
        if dfct is not None:
            # read the defect frame
            dfct = defect.MccdDefect.read(dfct)

        # Get windows from hdriver
        setup = cl.get_value("setup", "display current hdriver window settings", False)

        if setup:
            drurl = cl.get_value(
                "drurl", "URL for driver windows", "http://192.168.1.2:5100"
            )
        else:
            drurl = None

        # Some settings for the image plots
        cmap = cl.get_value(
            "cmap", "colour map to use ['none' for mpl default]", "Greys"
        )
        cmap = None if cmap == "none" else cmap

        imwidth = cl.get_value(
            "imwidth", "image plot width [inches, 0 for default]", 0.0, 0.0
        )

        imheight = cl.get_value(
            "imheight", "image plot height [inches, 0 for default]", 0.0, 0.0
        )

        memory = cl.get_value(
            "memory", "remember finaly plot dimensions for next time?", True
        )

        # define the display intensities
        msub = cl.get_value("msub", "subtract median from each window?", True)

        iset = cl.get_value(
            "iset",
            "set intensity a(utomatically)," " d(irectly) or with p(ercentiles)?",
            "a",
            lvals=["a", "d", "p"],
        )
        iset = iset.lower()

        plo, phi = 5, 95
        ilo, ihi = 0, 1000
        if iset == "d":
            ilo = cl.get_value("ilo", "lower intensity limit", 0.0)
            ihi = cl.get_value("ihi", "upper intensity limit", 1000.0)
        elif iset == "p":
            plo = cl.get_value(
                "plo", "lower intensity limit percentile", 5.0, 0.0, 100.0
            )
            phi = cl.get_value(
                "phi", "upper intensity limit percentile", 95.0, 0.0, 100.0
            )

        # region to plot
        for i, cnam in enumerate(ccds):
            nxtot, nytot, nxpad, nypad = ccdinf[cnam]
            if i == 0:
                xmin, xmax = float(-nxpad), float(nxtot + nxpad + 1)
                ymin, ymax = float(-nypad), float(nytot + nypad + 1)
            else:
                xmin = min(xmin, float(-nxpad))
                xmax = max(xmax, float(nxtot + nxpad + 1))
                ymin = min(ymin, float(-nypad))
                ymax = max(ymax, float(nytot + nypad + 1))

        xlo = cl.get_value("xlo", "left-hand X value", xmin, xmin, xmax, enforce=False)
        xhi = cl.get_value("xhi", "right-hand X value", xmax, xmin, xmax, enforce=False)
        ylo = cl.get_value("ylo", "lower Y value", ymin, ymin, ymax, enforce=False)
        yhi = cl.get_value("yhi", "upper Y value", ymax, ymin, ymax, enforce=False)

        # many parameters for profile fits, although most are not
        # plotted by default
        profit = cl.get_value("profit", "do you want accumulating profile fits?", False)

        #            nxf = cl.get_value("nxf", "maximum number of fit panels in X", 3)
        method = cl.get_value(
            "method", "fit method g(aussian) or m(offat)", "m", lvals=["g", "m"]
        )
        if method == "m":
            beta = cl.get_value(
                "beta", "initial exponent for Moffat fits", 5.0, 0.5, 20.0
            )
        else:
            beta = 0.0
        fwhm_min = cl.get_value(
            "fwhm_min", "minimum FWHM to allow [unbinned pixels]", 1.5, 0.01
        )
        fwhm = cl.get_value(
            "fwhm",
            "initial FWHM [unbinned pixels] for profile fits",
            6.0,
            fwhm_min,
        )
        shbox = cl.get_value(
            "shbox",
            "half width of box for initial location" " of target [unbinned pixels]",
            11.0,
            2.0,
        )
        smooth = cl.get_value(
            "smooth",
            "FWHM for smoothing for initial object" " detection [binned pixels]",
            6.0,
        )
        fhbox = cl.get_value(
            "fhbox",
            "half width of box for profile fit" " [unbinned pixels]",
            21.0,
            3.0,
        )
        hmin = cl.get_value("hmin", "minimum peak height to accept the fit", 50.0)
        if bias is None or msub:
            cl.set_default("read", -1.0)
        read = cl.get_value("read", "readout noise, RMS ADU", -1.0)
        gain = cl.get_value("gain", "gain, ADU/e-", 1.0)
        thresh = cl.get_value("thresh", "number of RMS to reject at", 4.0)

        fwnmax = cl.get_value(
            "fwnmax", "maximum number of frame to buffer FWHM [0 to ignore]", 100, 0
        )

        fwymax = cl.get_value(
            "fwymax",
            "maximum FWHM value to plot on Y-axis [unbinned pixels]",
            10.0,
            2.0,
        )

        fwwidth = cl.get_value(
            "fwwidth", "FWHM plot width [inches, 0 for default]", 0.0, 0.0
        )

        fwheight = cl.get_value(
            "fwheight", "FWHM plot height [inches, 0 for default]", 0.0, 0.0
        )

        cl.set_default("fcosmic", False)
        fcosmic = cl.get_value("fcosmic", "find cosmic rays", False)
        # disable if lacosmic is not available
        if lacosmic is None and fcosmic:
            warnings.warn("lacosmic not available, disabling cosmic ray removal")
            fcosmic = False

        if fcosmic:
            crthresh = cl.get_value(
                "crthresh", "threshold S-to-N for cosmic ray detection", 5.0, 0.0
            )
            nbthresh = cl.get_value(
                "nbthresh", "threshold S-to-N for neighbour pixels", 3.0, 0.0
            )
            sathresh = cl.get_value("sathresh", "saturation threshold", 65535.0)
            readout = cl.get_value("readout", "readout noise, RMS ADU", 2.5, 0.0)

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

    # Phew. We finally have all the inputs and now can now display stuff.

    # track which CCDs have been plotted at least once for the profile fits
    plotted = np.array(len(ccds) * [False])
    not_selected = True

    # img_accum and (later) fit_accum and targ_accum are long-lived
    # objects passed and stored by reference to plot updating classes
    img_accum = len(ccds) * [None]
    targ_accum = None
    first_image = True
    first_fwhm = True

    # Now onto the animated plots.
    with spooler.data_source(source, resource, first, full=False) as spool:

        nframe, total_time = 0, 0.0
        for mccd in spool:

            if server_or_local:
                # Handle the waiting game ... some awkward stuff
                # involving updating on a cycle faster than twait to
                # make the plots more responsive, if twait is long.
                if first_image and first_fwhm:
                    updaters = ()
                elif first_fwhm:
                    updaters = (imanager.update,)
                else:
                    updaters = (imanager.update, fwhmmanager.update)

                give_up, try_again, total_time = spooler.hang_about(
                    mccd, twait, tmax, total_time, updaters, tupdate
                )

                if give_up:
                    print("nrtplot stopped")
                    break
                elif try_again:
                    continue

            # Trim the frames: ULTRACAM windowed data has bad columns
            # and rows on the sides of windows closest to the readout
            # which can badly affect reduction. This option strips
            # them.
            if trim:
                hcam.ccd.trim_ultracam(mccd, ncol, nrow)

            # indicate progress
            tstamp = Time(mccd.head["TIMSTAMP"], format="isot", precision=3)
            nfrm = mccd.head.get("NFRAME", nframe + 1)
            print(
                f'{nfrm}, utc= {tstamp.iso} ({"ok" if mccd.head.get("GOODTIME", True) else "nok"}), ',
                end="",
            )

            # accumulate errors
            emessages = []

            # bias level checks
            if lowlevel != 0.0:
                median = mccd.get_num(0).get_num(0).median()
                if median < lowlevel:
                    emessages.append(
                        "** low bias level, median vs limit: {:.1f} vs {:.1f}".format(
                            median, lowlevel
                        )
                    )

            if highlevel != 0.0:
                try:
                    median = mccd.get_num(0).get_num(1).median()
                except:
                    median = mccd.get_num(0).get_num(0).median()

                if median > highlevel:
                    emessages.append(
                        "** high bias level, median vs limit: {:.1f} vs {:.1f}".format(
                            median, lowlevel
                        )
                    )

            if nframe == 0:

                # get the bias, dark, flat, fringe map and fringe pair files
                # into shape first time through

                if bias is not None:
                    # crop the bias on the first frame only
                    bias = bias.crop(mccd)
                    bexpose = bias.head.get("EXPTIME", 0.0)
                else:
                    bexpose = 0.0

                if dark is not None:
                    # crop the dark on the first frame only
                    dark = dark.crop(mccd)

                if flat is not None:
                    # crop the flat on the first frame only
                    flat = flat.crop(mccd)

                if fmap is not None:
                    # crop the fringe map and pair file
                    fmap = fmap.crop(mccd)
                    fpair = fpair.crop(mccd, nhalf)

            if setup:
                # Get setup windows from hdriver / udriver. Fair bit of
                # error checking needed. 'got_windows' indicates if
                # anything useful found, 'swindows' is a list of
                # (llx,lly,nx,ny) tuples if somthing is found.
                try:
                    r = requests.get(drurl, timeout=0.2)

                    if r.text.strip() == "No valid data available":
                        emessages.append(
                            f"** bad return from hdriver = {r.text.strip()}"
                        )
                        got_windows = False

                    elif r.text.strip() == "fullframe":
                        # to help Stu out a bit, effectively just
                        # ignore this one
                        got_windows = False

                    else:
                        # OK, got something
                        got_windows = True
                        lines = r.text.split("\r\n")
                        xbinh, ybinh, nwinh = lines[0].split()
                        xbinh, ybinh, nwinh = int(xbinh), int(ybinh), int(nwinh)
                        swindows = []
                        for line in lines[1 : nwinh + 1]:
                            llxh, llyh, nxh, nyh = line.split()
                            swindows.append((int(llxh), int(llyh), int(nxh), int(nyh)))

                        if nwinh != len(swindows):
                            emessages.append(
                                f"** expected {nwinh} windows from"
                                " (u/h)driver but got {len(swindows)}"
                            )
                            got_windows = False

                except (
                    requests.exceptions.ConnectionError,
                    socket.timeout,
                    requests.exceptions.Timeout,
                ) as err:
                    emessages.append(" ** hdriver error: {!r}".format(err))
                    got_windows = False

            else:
                got_windows = False

            # wind through the CCDs to display, accumulating stuff
            # to send to the plot manager
            message = ""

            skipped = True
            for nc, cnam in enumerate(ccds):
                ccd = mccd[cnam]

                if plotall or ccd.is_data():
                    # "is_data" indicates genuine data as opposed to junk
                    # that results from nskip > 0.
                    plotted[nc] = True

                    if fcosmic:
                        # save a copy for later in order to mask
                        # saturated stars.
                        raw = ccd.copy()

                    # subtract the bias
                    if bias is not None:
                        ccd -= bias[cnam]

                    if dark is not None:
                        dexpose = dark.head["EXPTIME"]
                        cexpose = ccd.head["EXPTIME"]
                        scale = (cexpose - bexpose) / dexpose
                        ccd -= scale * dark[cnam]

                    # divide out the flat
                    if flat is not None:
                        ccd /= flat[cnam]

                    # Remove fringes
                    if fmap is not None and cnam in fmap and cnam in fpair:
                        fscale = fpair[cnam].scale(ccd, fmap[cnam], nhalf, rmin, rmax)
                        ccd -= fscale * fmap[cnam]

                    # cosmic ray detection
                    if fcosmic:
                        ncosmic = 0
                        for wind, rwind in zip(ccd.values(), raw.valued()):
                            data = wind.data
                            smask = rwind.data >= sathresh
                            cleaned, crmask = lacosmic.lacosmic(
                                wind.data,
                                crthresh,
                                nbthresh,
                                mask=smask,
                                effective_gain=gain,
                                readnoise=readout * gain,
                            )
                            ncosmic += len(crmask[crmask])

                    if msub:
                        # subtract median from each window
                        for wind in ccd.values():
                            wind -= wind.median()

                    # set intensity limits
                    if iset == "p":
                        # Set intensities from percentiles
                        if xlo is not None and xhi is not None:
                            pxlo = min(xlo, xhi)
                            pxhi = max(xlo, xhi)
                        else:
                            pxlo, pxhi = None, None
                        if ylo is not None and yhi is not None:
                            pylo = min(ylo, yhi)
                            pyhi = max(ylo, yhi)
                        else:
                            pylo, pyhi = None, None
                        vmin, vmax = ccd.percentile((plo, phi), pxlo, pxhi, pylo, pyhi)

                    elif iset == "a":
                        # Set intensities from min/max range
                        vmin, vmax = ccd.min(), ccd.max()

                    elif iset == "d":
                        vmin, vmax = ilo, ihi

                    content = [ccd, vmin, vmax]

                    if got_windows:
                        content.append(swindows)
                    else:
                        content.append(None)

                    if dfct is not None and cnam in dfct:
                        content.append(dfct[cnam])
                    else:
                        content.append(None)

                    # Save the content to send to the plot updater
                    img_accum[nc] = content

                    # accumulate string of image scalings
                    if nc:
                        message += ", "
                    message += f"ccd {cnam}: {vmin:.1f}, {vmax:.1f}, exp: {mccd.head['EXPTIME']:.4f}"
                    if fcosmic:
                        message += f", nc: {ncosmic}"
                    skipped = False

            # Print messages
            if skipped:
                print(f"{message}skipped")
            else:
                print(f"{message}")
            for emessage in emessages:
                print(emessage)

            # at this point "img_accum" contains a list of lists, each of
            # which consists of:
            #
            # [ccd, vmin, vmax, swindows, dfct]
            #
            # i.e. the CCD, intensity range, setup windows and defects, or
            # "None" if the CCD has yet to be initialised due to nskip

            if profit and not_selected and plotted.all():
                # Finally have at least one proper exposure of all CCDs
                # from which stars can be selected.
                not_selected = False

                # Ready to make first plot of all CCDs
                img_fig, img_axs = setup_images(
                    len(ccds),
                    nx,
                    "Profile fit selection",
                    xlo,
                    xhi,
                    ylo,
                    yhi,
                    imwidth,
                    imheight,
                )
                for ax, cnam, content in zip(img_axs, ccds, img_accum):
                    disp_ccd(ax, cnam, cmap, content, False)

                # Cursor selection routine. Comes back with a set of targets
                cselect = ProfitCursorSelect(
                    img_fig,
                    img_axs,
                    ccds,
                    img_accum,
                    shbox,
                    fwhm,
                    beta,
                    method,
                    smooth,
                    fhbox,
                    hmin,
                    fwhm_min,
                    read,
                    gain,
                    thresh,
                )
                plt.show()

                targets = cselect.targets
                if len(targets) == 0:
                    # nothing selected, give up on
                    # profile fits.
                    profit = False

            if profit and not not_selected:
                # If profile fitting is being carried out, then this
                # section is implemented for each new image.

                # First update search box positions. In normal
                # operation this moves the boxes to match the target
                # positions previously fitted.
                for targ in targets:
                    targ.update_sbox()

                # Carry out fits
                pfits = prof_fit(
                    ccds,
                    targets,
                    img_accum,
                    method,
                    smooth,
                    fhbox,
                    hmin,
                    fwhm_min,
                    read,
                    gain,
                    thresh,
                )

                # if we have two (and only two) valid fits, then find the angle between them
                if len(pfits) == 2 and pfits[0] is not None and pfits[1] is not None:
                    dx = pfits[1][0].x - pfits[0][0].x
                    dy = pfits[1][0].y - pfits[0][0].y
                    theta = Angle(np.arctan2(dy, dx), unit=u.rad).to(u.deg)
                    theta = theta.wrap_at(360*u.deg).to_value()
                    print(f"  From target 1 to 2: \u03F4 = {theta:.2f}\u00B0 (anti-clockwise from horizontal)")

                if first_fwhm:
                    # set up the accumulator object now because we didn't know
                    # how many targets the user would choose at the start
                    fit_accum = [([], []) for n in range(len(pfits))]
                    targ_accum = len(pfits) * [None]

                # add in latest results
                for n, pfit in enumerate(pfits):
                    if pfit is not None:
                        fwhm = pfit[0].fwhm
                        fit_accum[n][0].append(nfrm)
                        fit_accum[n][1].append(fwhm)
                        targ_accum[n] = pfit[0]

            if not profit or not not_selected:

                # The image animation
                if first_image:
                    # one-off setup
                    first_image = False

                    # Create the image plot manager
                    imanager = ImageManager(
                        ccds,
                        nx,
                        imwidth,
                        imheight,
                        xlo,
                        xhi,
                        ylo,
                        yhi,
                        cmap,
                        img_accum,
                        targ_accum,
                        shbox,
                        fwhm,
                        beta,
                        method,
                        smooth,
                        fhbox,
                        hmin,
                        fwhm_min,
                        read,
                        gain,
                        thresh,
                    )

                    if memory:
                        # define object to handle adjusting the figure size
                        cleanup = CleanUp(
                            "HIPERCAM_ENV", ".hipercam", command, args, imanager.fig
                        )

                        # and register with signal
                        signal(SIGINT, cleanup)

                    plt.show(block=False)
                    plt.pause(0.1)

                else:
                    # update plot
                    imanager.update()

                if profit and fwnmax > 1:
                    # The fwhm animation
                    if first_fwhm:
                        # one-off setup
                        first_fwhm = False

                        # Create the FWHM plot manager
                        fwhmmanager = FwhmManager(
                            fit_accum,
                            nfrm - 1,
                            nfrm + fwnmax + 1,
                            fwymax,
                            fwwidth,
                            fwheight,
                        )

                        plt.show(block=False)
                        plt.pause(0.1)

                    else:
                        # update the fwhm plot manager
                        fwhmmanager.update()

            if not (profit and not_selected) and pause > 0.0:
                # pause between frames
                if pause > 2 * tupdate:
                    # run the updaters if pause is a little
                    # long. Designed to make the plot more responsive
                    nupdate = int(pause / tupdate)
                    for n in range(nupdate):
                        time.sleep(pause / nupdate)
                        if not first_image:
                            imanager.update()
                        if not first_fwhm:
                            fwhmmanager.update()
                else:
                    # just pause because the updates will come
                    time.sleep(pause)

            # update the frame number
            nframe += 1

    # At then end the plots will close automatically,
    # so block this by asking the user to confirm
    reply = input("\nHit carriage return to close the plot(s): ")

    if memory:
        cleanup.finish()


# From here is support code not visible outside


class CleanUp:
    """
    Handle adjusting plot width defaults in a way that can be used
    as a handler for signal
    """

    def __init__(self, hipenv, hipconf, command, args, fig):
        self.hipenv = hipenv
        self.hipconf = hipconf
        self.command = command
        self.args = args
        self.fig = fig

    def __call__(self, signal_recieved, frame):
        self.finish()
        sys.exit(1)

    def finish(self):
        with Cline(self.hipenv, self.hipconf, self.command, self.args) as cl:
            width = self.fig.get_figwidth()
            height = self.fig.get_figheight()
            cl.register("imwidth", Cline.LOCAL, Cline.HIDE)
            cl.register("imheight", Cline.LOCAL, Cline.HIDE)
            cl.set_default("imwidth", width)
            cl.set_default("imheight", height)
            print(
                f"\nImage display size updated to width, height = {width}, {height} [memory=True]"
            )


def setup_images(nccd, nx, title, xlo, xhi, ylo, yhi, width, height):
    """
    Sets up the figure and axes for the images display.

    Arguments::

       nccd : int
          number of CCDs to display

       nx : int
          number of panels in the X direction.

       title : str
          title of figure

       xlo, xhi, ylo, yhi : floats
          the unbinned pixel eange to display

       width : float
          width, inches

       height : float
          height, inches

    Returns (fig, axs)

    The Figure and Axes, one per CCD
    """

    # Create the figure for the image plot
    if width > 0 and height > 0:
        fig = plt.figure(title, figsize=(width, height))
    else:
        fig = plt.figure(title)

    # Define config of images
    ny = nccd // nx if nccd % nx == 0 else nccd // nx + 1

    # Create the axes
    for ind in range(nccd):
        if ind == 0:
            ax0 = ax = fig.add_subplot(ny, nx, ind + 1, aspect="equal")
            ax.set_xlim(xlo, xhi)
            ax.set_ylim(ylo, yhi)
            axs = [ax]
        else:
            ax = fig.add_subplot(
                ny, nx, ind + 1, sharex=ax0, sharey=ax0, aspect="equal"
            )
            axs.append(ax)
        ax.set_aspect("equal")
        for tick in ax.get_yticklabels():
            tick.set_rotation(90)

    return (fig, axs)


class ImageManager:
    """Class to control the image animation

    Basic idea is that initialising it stores all the very many inputs
    as attributes and creates the plot returning structures containing
    animated artists which get updated per frame. Allows re-scaling of
    intensity levels of the images, plotting of defects and setup
    windows and also the search box and fitted positions of any
    profile fit targets.

    """

    def __init__(
        self,
        cnams,
        nx,
        width,
        height,
        xlo,
        xhi,
        ylo,
        yhi,
        cmap,
        img_accum,
        targ_accum,
        shbox,
        fwhm,
        beta,
        method,
        smooth,
        fhbox,
        hmin,
        fwhm_min,
        read,
        gain,
        thresh,
    ):
        """Initialises the plot. Arguments:

          cnams : list of strings
             CCD names to plot

          nx : int
             max number of panels in X

          width : float
             plot width [inches]

          height : float
             plot height [inches]

          xlo : float
             left-hand limit of region to determine percentiles and for display

          xhi : float
             right-hand limit of region to display

          ylo : float
             lower limit of region to display

          yhi : float
             upper limit of region to display

          cmap : str  | None
             Matplotlib colour map. None gives default


          img_accum : list of lists
             one list per CCD being displayed. Contains the information
             required for the image plot. Per CCD the list should be:
             [ccd, vmin, vmax, swindows, dfcts] or None if the CCD is
             being skipped. swindows -- setup windows -- can be None.
             dfcts -- defects -- can also be None. A reference to this
             is stored internally so that changes to it outside the object
             are reflected internally. This allows update to work without
             arguments.

          targ_accum : list of Fpars
             one tuple per target being profile fitted, containg all
             info needed to make plots of the fit etc. This can be
             None if no fit was made. See Fpar.fit for what the tuples
             contain.  A reference to this is stored internally so
             that changes to it outside the object are reflected
             internally. This allows update to work without arguments.

        shbox, fwhm, beta, method, smooth, fhbox, hmin, fwhm_min, read,
        gain, thresh are profile fitting parameters described elsewhere.
        They are sent to a cursor selection routine.

        """

        # basic check
        assert len(img_accum) == len(cnams)

        # create figure and axes
        self.fig, self.axs = setup_images(
            len(cnams), nx, "CCD display", xlo, xhi, ylo, yhi, width, height
        )

        self.cnv = self.fig.canvas
        self._bg = None
        self.img_accum = img_accum
        self.targ_accum = targ_accum

        # nrtplot inputs
        self.cnams = cnams
        self.xlo = xlo
        self.xhi = xhi
        self.ylo = ylo
        self.yhi = yhi
        self.cmap = cmap

        # list of dictionaries of artists for each image panel
        self.img_artists = []

        # now we actually create the artists
        for ax, cnam, content in zip(self.axs, self.cnams, self.img_accum):
            if content is None:
                # do nothing much other than add a placeholder
                self.img_artists.append(None)
            else:
                # plot the CCD, storing references to the animated artists
                self.img_artists.append(self._disp_ccd(ax, cnam, content))

        # now for the fit targets
        self.targ_artists = self._disp_targs()

        # grab the background on every draw
        self.cid = self.cnv.mpl_connect("draw_event", self.on_draw)

        # profile fitter
        self.cselect = OntheflyCursorSelect(
            self.fig,
            self.axs,
            self.cnams,
            img_accum,
            shbox,
            fwhm,
            beta,
            method,
            smooth,
            fhbox,
            hmin,
            fwhm_min,
            read,
            gain,
            thresh,
        )

    def on_draw(self, event):
        """Callback to register with 'draw_event'."""
        cnv = self.cnv
        if event is not None and event.canvas != cnv:
            raise RuntimeError
        self._bg = cnv.copy_from_bbox(cnv.figure.bbox)
        self._draw_animated()

    def _draw_animated(self):
        """Draws all of the animated artists for the image plot self.artists
        is a list of dictionaries of lists ...

        """

        # need to know some of the details behind the dictionaries
        # returned by mpl
        for dartist in self.img_artists:

            # dartist is a dictionary as returned by disp_ccd

            if dartist is not None:

                # 'ccd' has values which are all lists of artists
                for lartist in dartist["ccd"].values():
                    for artist in lartist:
                        self.fig.draw_artist(artist)

                # so does 'defects'
                for lartist in dartist["defects"].values():
                    for artist in lartist:
                        self.fig.draw_artist(artist)

                # but 'swindows' is a plain list of artists
                for artist in dartist["swindows"]:
                    self.fig.draw_artist(artist)

        # targets
        for artist in self.targ_artists:
            if artist is not None:
                self.fig.draw_artist(artist)

    def update(self):
        """Updating routine."""
        # update the frames used by the cursor picker
        self.cselect.update(self.img_accum)
        # now update / create the artists

        # first the images
        lartists = []
        for ax, cnam, content, artists in zip(
            self.axs, self.cnams, self.img_accum, self.img_artists
        ):
            if content is None:
                # just pass old artists through
                lartists.append(artists)
            else:
                # plot the CCD, return with the updated animated artists
                lartists.append(self._disp_ccd(ax, cnam, content, artists))
        self.img_artists = lartists

        # then the profile fit targets
        self.targ_artists = self._disp_targs(self.targ_artists)

        cnv = self.cnv
        fig = self.fig
        if self._bg is None:
            self.on_draw(None)
        else:
            # restore the background
            cnv.restore_region(self._bg)
            self._draw_animated()
            cnv.blit(fig.bbox)
        cnv.flush_events()

    def _disp_ccd(self, ax, cnam, content, artists=None):
        """Displays a CCD ccd, name cnam, in Axes ax.

        If "artists" is None, a dictionary will be created (see
        disp_ccd for details). Otherwise it is assumed to be such a
        dcitionary resulting from an earlier run and will be updated.

        """

        if artists is None:

            # in this case we are setting up for the first time
            artists = disp_ccd(ax, cnam, self.cmap, content, True)

            # need to re-draw to avoid irritating distorted image
            self.cnv.draw()

        else:

            # unpack the new content
            ccd, vmin, vmax, swindows, dfct = content

            # this is the "usual" post-setup case where we just update
            # the artists

            # the CCD window data
            mpl.pCcd(
                ax,
                ccd,
                "d",
                dlo=vmin,
                dhi=vmax,
                tlabel=f"CCD {cnam}",
                cmap=self.cmap,
                animated=True,
                artists=artists["ccd"],
            )

            # the setup windows. Can change in number as well as size
            # so delete old ones (if there are any) and re-create each
            # time remove old setup windows
            swins = artists["swindows"]
            for artist in swins:
                artist.remove()
            artists["swindows"] = []

            if swindows is not None:

                # draw new ones, keeping refs to objects for next time
                swins = []
                for llxh, llyh, nxh, nyh in swindows:
                    (box,) = ax.plot(
                        [
                            llxh - 0.5,
                            llxh + nxh - 0.5,
                            llxh + nxh - 0.5,
                            llxh - 0.5,
                            llxh - 0.5,
                        ],
                        [
                            llyh - 0.5,
                            llyh - 0.5,
                            llyh + nyh - 0.5,
                            llyh + nyh - 0.5,
                            llyh - 0.5,
                        ],
                        "--",
                        color=COL_SETUP,
                        animated=True,
                    )
                    swins.append(box)
                artists["swindows"] = swins

        return artists

    def _disp_targs(self, artists=None):
        """Displays target as (a) a search box and (b) a circle of radius
        equal to the FWHM and a dot for the target.

        If "artists" is None, a list will be created and
        returned containing the animated artists.
        Otherwise it is assumed to be such a list resulting from an
        earlier run and will be updated.

        """

        if self.targ_accum is None:
            return []

        # mapping from CCD name to axes
        ccd2axes = dict(zip(self.cnams, self.axs))

        if artists is None:

            artists = []
            for ntarg, fpar in enumerate(self.targ_accum):
                if fpar is not None:
                    ax = ccd2axes[fpar.cnam]
                    artists += fpar.plot(ax, True)
                else:
                    artists += [None, None, None]

            # re-draw
            self.cnv.draw()

        else:

            for ntarg, fpar in enumerate(self.targ_accum):
                if fpar is not None:
                    if artists[3 * ntarg] is None:
                        # The artists have yet to be made
                        ax = ccd2axes[fpar.cnam]
                        artists[3 * ntarg : 3 * ntarg + 3] = fpar.plot(ax, True)
                    else:
                        # they exist already and need updating
                        dot, circ, sbox = artists[3 * ntarg : 3 * ntarg + 3]
                        dot.set_data(fpar.x, fpar.y)
                        circ.set_center((fpar.x, fpar.y))
                        circ.set_radius(fpar.fwhm)
                        xlo, xhi, ylo, yhi = fpar.region()
                        sbox.set_data(
                            [xlo, xhi, xhi, xlo, xlo], [ylo, ylo, yhi, yhi, ylo]
                        )

        return artists


def disp_ccd(ax, cnam, cmap, content, animated):
    """Displays a CCD ccd, name cnam, in Axes ax.

    If used for movie-style display, set animated=True.

    Returns dictionary of artists with keys:

      'ccd' : dictionary as returned by pCcd
      'defects' : dictionary as returned by pCcdDefect. Can be just {}
      'swindows' : list of Artists, can be just []

    """

    # unpack the new content
    ccd, vmin, vmax, swindows, dfct = content

    artists = {}

    # plot the ccd
    artists["ccd"] = mpl.pCcd(
        ax,
        ccd,
        "d",
        dlo=vmin,
        dhi=vmax,
        tlabel=f"CCD {cnam}",
        cmap=cmap,
        animated=animated,
    )[2]

    artists["defects"] = mpl.pCcdDefect(ax, dfct, animated)

    # plot setup windows
    if swindows is None:
        swins = []
    else:
        swins = []
        for llxh, llyh, nxh, nyh in swindows:
            (box,) = ax.plot(
                [
                    llxh - 0.5,
                    llxh + nxh - 0.5,
                    llxh + nxh - 0.5,
                    llxh - 0.5,
                    llxh - 0.5,
                ],
                [
                    llyh - 0.5,
                    llyh - 0.5,
                    llyh + nyh - 0.5,
                    llyh + nyh - 0.5,
                    llyh - 0.5,
                ],
                "--",
                color=COL_SETUP,
                animated=animated,
            )
            swins.append(box)
    artists["swindows"] = swins
    return artists


class FwhmManager:
    """Class to start and control the FWHM history animation"""

    def __init__(self, accum, xlo, xhi, fwmax, width, height):
        """
        Initialises the FWHM history plot. Arguments:

          accum : list of of tuples
             the lists consists of a sequence of tuples with the frame number
             and results for all targets (nframe, pfits), with pfits a seies
             of fpar objects (can be None if no fit made)

          xlo : float
             left-hand X -limit (inital) for plot

          xhi : float
             right-hand X -limit (inital) for plot

          fwmax : float
             maximum height for FWHM plot in pixels

          width : float
             plot width [inches]

          height : float
             plot height [inches]

        """

        self.accum = accum

        # Create the figure
        if width > 0 and height > 0:
            self.fig = plt.figure("FWHM vs frame number", figsize=(width, height))
        else:
            self.fig = plt.figure("FWHM vs frame number")

        self.cnv = self.fig.canvas
        self.ax = self.fig.add_subplot(animated=True)
        self._bg = None

        # nrtplot inputs
        self.fwmax = fwmax

        # derived
        self.xlo = xlo
        self.xhi = xhi
        self.ylo = 0
        self.yhi = fwmax

        self.ax.set_xlim(self.xlo, self.xhi)
        self.ax.set_ylim(self.ylo, self.yhi)
        self.ax.set_xlabel("Frame number")
        self.ax.set_ylabel("FWHM [unbinned pixels]")

        # Now plot and create artists for each FWHM plot
        self.fwhm_artists = []
        for ntarg, (xs, ys) in enumerate(self.accum):
            (fwhm,) = self.ax.plot(xs, ys, animated=True, label=f"{ntarg+1}")
            self.fwhm_artists.append(fwhm)
        self.ax.legend()
        #        self.cnv.draw()

        # grab the background on every draw
        self.cid = self.cnv.mpl_connect("draw_event", self.on_draw)

    def on_draw(self, event):
        """Callback to register with 'draw_event'."""
        cnv = self.cnv
        if event is not None and event.canvas != cnv:
            raise RuntimeError
        self._bg = cnv.copy_from_bbox(cnv.figure.bbox)
        self._draw_animated()

    def _draw_animated(self):
        # draw all of the animated artists. NB draw
        # Axes before the lines if you want to see anything!
        self.fig.draw_artist(self.ax)

        for artist in self.fwhm_artists:
            if artist is not None:
                self.fig.draw_artist(artist)

    def update(self):
        """updating routine. Uses current value of self.accum which
        should be a list of Fpar objects for each target.

        """

        xmax = None
        for ntarg, (xs, ys) in enumerate(self.accum):
            self.fwhm_artists[ntarg].set_data(xs, ys)
            if len(xs):
                xmax = max(xs) if xmax is None else max(xmax, max(xs))

        if xmax > self.xhi - 1:
            # update axes limits
            xadd = xmax - self.xhi + 1
            self.xlo += xadd
            self.xhi += xadd
            self.ax.set_xlim(self.xlo, self.xhi)

        cnv = self.cnv
        fig = self.fig
        if self._bg is None:
            self.on_draw(None)
        else:
            # restore the background
            cnv.restore_region(self._bg)
            self._draw_animated()
            cnv.blit(fig.bbox)
        cnv.flush_events()


class ProfitCursorSelect:
    """
    Selection of targets to profile fit for animation
    """

    def __init__(
        self,
        fig,
        axs,
        ccds,
        img_accum,
        shbox,
        fwhm,
        beta,
        method,
        smooth,
        fhbox,
        hmin,
        fwhm_min,
        read,
        gain,
        thresh,
    ):
        """
        fig : the Figure
        axs : the Axes, one per CCD
        img_accum : image content data
        shbox : falf width search box
        """

        self.cnv = fig.canvas
        self.axs = axs
        self.ccds = ccds
        self.img_accum = img_accum
        self.shbox = shbox
        self.fwhm = fwhm
        self.beta = beta
        self.method = method
        self.smooth = smooth
        self.fhbox = fhbox
        self.hmin = hmin
        self.fwhm_min = fwhm_min
        self.read = read
        self.gain = gain
        self.thresh = thresh

        # accumulating target list
        self.targets = []
        self.ntarg = 0
        print(
            "\nUse the SpaceBar to select targets for profile fitting; q(uit) to exit selection and start animation"
        )
        self.kp = self.cnv.mpl_connect("key_press_event", self._keyPressEvent)

    def _keyPressEvent(self, event):
        """
        Where stuff is done
        """

        if event.key == " " and event.inaxes is not None:

            for ax, cnam, content in zip(self.axs, self.ccds, self.img_accum):

                # search to identify the CCD
                if event.inaxes == ax:

                    # clicked inside an Axes
                    x, y = event.xdata, event.ydata

                    # extract the CCD
                    ccd = content[0]

                    # check that the position is inside a window
                    wnam = ccd.inside(x, y, 2)

                    if wnam is not None:
                        # store the position, Window label, target
                        # number, box size fwhm, beta
                        fpar = Fpar(cnam, wnam, x, y, self.shbox, self.fwhm, self.beta)
                        results, message = fpar.fit(
                            ccd,
                            self.method,
                            self.smooth,
                            self.fhbox,
                            self.hmin,
                            self.fwhm_min,
                            self.read,
                            self.gain,
                            self.thresh,
                        )
                        if results is not None:
                            # fitted OK
                            self.ntarg += 1
                            self.targets.append(fpar)
                            print(
                                f"\n   target {self.ntarg} added at initial x,y = "
                                f"{x:.2f}, {y:.2f} in CCD {cnam}, window {wnam}, refined by fit with results:"
                            )
                            print(f"   {message}")
                            fpar.plot(ax)
                            self.cnv.draw()
                        else:
                            print(
                                f"\n   ** fit failed for position x,y = {x:.2f}, {y:.2f} in CCD {cnam}, window {wnam}; no new target added"
                            )
                            print(f"   ** fit message = {message}")


class OntheflyCursorSelect:
    """
    Selection of targets for on-the-fly fits
    """

    def __init__(
        self,
        fig,
        axs,
        ccds,
        img_accum,
        shbox,
        fwhm,
        beta,
        method,
        smooth,
        fhbox,
        hmin,
        fwhm_min,
        read,
        gain,
        thresh,
    ):
        """
        fig : the Figure
        axs : the Axes, one per CCD
        img_accum : image content data
        shbox : falf width search box
        """

        self.cnv = fig.canvas
        self.axs = axs
        self.ccds = ccds
        self.img_accum_old = None
        # have to copy here to avoid reference-passing auto-updating
        # the images because we have to work one frame behind
        self.img_accum = img_accum.copy()
        self.shbox = shbox
        self.fwhm = fwhm
        self.beta = beta
        self.method = method
        self.smooth = smooth
        self.fhbox = fhbox
        self.hmin = hmin
        self.fwhm_min = fwhm_min
        self.read = read
        self.gain = gain
        self.thresh = thresh

        print("\n(use the SpaceBar to select targets for one-off profile fits)\n")
        self.kp = self.cnv.mpl_connect("key_press_event", self._keyPressEvent)

    def update(self, img_accum):
        self.img_accum_old = self.img_accum
        self.img_accum = img_accum.copy()

    def _keyPressEvent(self, event):
        """
        Where stuff is done
        """

        if (
            event.key == " "
            and event.inaxes is not None
            and self.img_accum_old is not None
        ):

            for ax, cnam, content in zip(self.axs, self.ccds, self.img_accum_old):

                if event.inaxes == ax:

                    # clicked inside an Axes
                    x, y = event.xdata, event.ydata

                    # extract the CCD
                    ccd = content[0]

                    # check that the position is inside a window
                    wnam = ccd.inside(x, y, 2)

                    if wnam is not None:
                        # store the position, Window label, target number,
                        # box size fwhm, beta
                        fpar = Fpar(cnam, wnam, x, y, self.shbox, self.fwhm, self.beta)
                        results, message = fpar.fit(
                            ccd,
                            self.method,
                            self.smooth,
                            self.fhbox,
                            self.hmin,
                            self.fwhm_min,
                            self.read,
                            self.gain,
                            self.thresh,
                        )
                        if results is not None:
                            # fitted OK
                            print(
                                f"   profile fit with initial x,y = "
                                f"{x:.2f}, {y:.2f} in CCD {cnam}, window {wnam}, [applies to previous frame]:"
                            )
                            print(f"   {message}\n")
                        else:
                            print(
                                f"\n   ** fit failed at position x,y = {x:.2f}, {y:.2f} in CCD {cnam}, window {wnam}"
                            )
                            print(f"   ** fit message = {message}\n")

    def _onpress(self, event):
        self.press = True

    def _onmove(self, event):
        if self.press:
            self.move = True

    def _onrelease(self, event):
        # only call onclick in special circumstances
        if self.press and not self.move:
            self._onclick(event)
        self.press = False
        self.move = False


class Fpar:

    """Class for representing a target for profile fits. Stores
    the CCD name and window name where it is located."""

    def __init__(self, cnam, wnam, x, y, shbox, fwhm, beta):
        """
        x,y used as initial target position and initial centre of search box
        """
        self.cnam = cnam
        self.wnam = wnam
        self.shbox = shbox
        self.x_cbox = self.x = x
        self.y_cbox = self.y = y
        self.fwhm = fwhm
        self.beta = beta

    def region(self):
        """
        Returns current search region
        """
        return (
            self.x_cbox - self.shbox,
            self.x_cbox + self.shbox,
            self.y_cbox - self.shbox,
            self.y_cbox + self.shbox,
        )

    def plot(self, ax, animated=False):
        """Plots the Fpar as a search box plus and a circle of radius FWHM and
        a dot for target. Returns with a list of three artists. The
        artists can be animated if they are be updated.
        """

        xlo, xhi, ylo, yhi = self.region()

        (sbox,) = ax.plot(
            [xlo, xhi, xhi, xlo, xlo],
            [ylo, ylo, yhi, yhi, ylo],
            color="b",
            animated=animated,
        )
        circ = Circle(
            (self.x, self.y), self.fwhm, fill=False, color="g", animated=animated
        )
        ax.add_patch(circ)
        (dot,) = ax.plot(self.x, self.y, ".g", animated=animated)
        return (dot, circ, sbox)

    def fit(self, ccd, method, smooth, fhbox, hmin, fwhm_min, read, gain, thresh):
        """Carries out a fit when passed a CCD. If successful,
        x, y, fwhm and beta will be altered but the centre
        of the search region will be unchanged.

        Returns (results, message) where results = (self, fradii, fdata, fok, vmin, vmax, r, f)
        if successful, None if not. "message" is a string summarising the result.

        """

        if read < 0:
            # Set a fairly high readout noise to start with in this case
            # to equalise weights a bit as we are going to cycle it with a
            # second fit
            read1 = 40.0
        else:
            read1 = read

        try:
            # extract search box data from the CCD.
            xlo, xhi, ylo, yhi = self.region()
            swind = ccd[self.wnam].window(xlo, xhi, ylo, yhi)

            # carry out initial search
            x, y, peak = swind.search(smooth, self.x, self.y, hmin, False)

            # now for a more refined fit. First extract region centred
            # on new position
            fwind = ccd[self.wnam].window(x - fhbox, x + fhbox, y - fhbox, y + fhbox)

            # crude estimate of sky background
            sky = np.percentile(fwind.data, 50)

            # uncertainties
            sigma = np.sqrt(read1**2 + np.maximum(0, fwind.data) / gain)

            # refine the Aperture position by fitting the profile
            (
                (sky, height, x, y, fwhm, beta),
                epars,
                (wfit, X, Y, chisq, nok, nrej, npar, nfev, message),
            ) = hcam.fitting.combFit(
                fwind,
                sigma,
                method,
                sky,
                peak - sky,
                x,
                y,
                self.fwhm,
                fwhm_min,
                False,
                self.beta,
                20.0,
                False,
                thresh,
            )

            if read < 0:
                # Estimate the read noise post fit and have another go ...
                sfac = np.sqrt(chisq / nok)
                resid = np.abs((fwind.data - wfit.data) / sigma)
                ok = (sigma > 0) & (resid < 4 * sfac)
                read2 = np.std((fwind.data - wfit.data)[ok])
                ok = (
                    (sigma > 0)
                    & ((fwind.data - sky) < 4 * read2)
                    & (np.abs(wfit.data - sky) < read2)
                )
                if len(sigma[ok]) > 5:
                    var = np.var(fwind.data[ok])
                    if var < sky / gain:
                        # variance less than expected
                        corr = gain * (var - sky)
                        fwind -= corr
                        extra = (
                            f", set read = 0 and applied bias correction of {corr:.1f}."
                        )
                        sky -= corr
                        read2 = 0
                    else:
                        # compute read noise to give correct variance in sky
                        read2 = np.sqrt(var - sky / gain)
                        extra = f", read = {read2:.2f}"

                    # uncertainties
                    sigma = np.sqrt(read2**2 + np.maximum(0, fwind.data) / gain)

                    # refine the Aperture position by fitting the profile
                    (
                        (sky, height, x, y, fwhm, beta),
                        epars,
                        (wfit, X, Y, chisq, nok, nrej, npar, nfev, message),
                    ) = hcam.fitting.combFit(
                        fwind,
                        sigma,
                        method,
                        sky,
                        (wfit.data - sky).max(),
                        x,
                        y,
                        fwhm,
                        fwhm_min,
                        False,
                        beta,
                        20.0,
                        False,
                        thresh,
                    )
                    message += extra
                else:
                    message += f', failed to measure "read"; used value = {read1}'

            if peak > hmin and ccd[self.wnam].distance(x, y) > 1:
                # update some initial parameters for next time
                if method == "g":
                    self.x, self.y, self.fwhm = x, y, fwhm
                elif method == "m":
                    self.x, self.y, self.fwhm, self.beta = x, y, fwhm, beta

                # values vs radial distance, plot range
                ok = sigma > 0
                R = np.sqrt((X - x) ** 2 + (Y - y) ** 2)
                vmin = min(sky, sky + height, fwind.min())
                vmax = max(sky, sky + height, fwind.max())
                extent = vmax - vmin
                vmin -= 0.05 * extent
                vmax += 0.05 * extent

                # line fit
                r = np.linspace(0, R.max(), 400)
                if method == "g":
                    alpha = 4 * np.log(2.0) / fwhm**2
                    f = sky + height * np.exp(-alpha * r**2)
                elif method == "m":
                    alpha = 4 * (2 ** (1 / beta) - 1) / fwhm**2
                    f = sky + height / (1 + alpha * r**2) ** beta

                # save the Fpar, radii, data, flags, plot range, fit x & y
                results = (self, R.flat, fwind.data.flat, ok.flat, vmin, vmax, r, f)

            else:
                results = None
                message += "\n*** below detection threshold. Parameters unchanged."

        except hcam.HipercamError as err:
            results = None
            message = f"fit failed: {err}"

        return (results, message)

    def update_sbox(self):
        """Moves centre of box to match target position"""
        self.x_cbox = self.x
        self.y_cbox = self.y


def prof_fit(
    ccds, targets, img_accum, method, smooth, fhbox, hmin, fwhm_min, read, gain, thresh
):
    """Carries out profile fitting over all targets. Returns container accumulating x,y lists for plotting
    the FWHM vs frame for each target.

    """

    # first a mapping from CCD name to the image content
    img_map = dict(zip(ccds, img_accum))

    # container for results
    pfits = len(targets) * [None]

    for ntarg, fpar in enumerate(targets):
        # carry out fits. Nothing happens if the CCD has not been
        # updated [content = None]
        content = img_map[fpar.cnam]
        if content is not None:
            ccd = content[0]
            results, message = fpar.fit(
                ccd, method, smooth, fhbox, hmin, fwhm_min, read, gain, thresh
            )
            print(f" target {ntarg+1}, {message}")
            pfits[ntarg] = results

    return pfits
