# Last edited on 2025-08-19 12:18:42 by stolfi

#
# PLOT FILES
#
# The program also writes for each {mask} a file "{mdir}/spectrum-plot.png"
# for plotting the mean spectra {savg[:im]} for all masks together.
# 
# The first two lines are "nb = {nb}" and "nm = {nm}".
# Then follow a line "masks = {masks[0]} {masks[1]} ... {masks[nm-1]}"
# with the names of the masks (e.g. "ch", "di", "bg").
# Then follows {nb} lines with format
#
#   "{bands[ib]} {illum[ib]} {wlen[ib]} {data[ib,0]} ... {data[ib,nm-1]}"
#
# where {data[ib,im]} is
#
#   "{savg[ib,im]} {smin[ib,im]} {smax[ib,im]} {tdev[ib,im]} {P[ib,0,im]}"
#
# and {tdev[ib,im]} is the magnitude of the residual not accounted for
# by the first component {P[ib,0,im]}, that is, {tdev[ib,im] =
# sqrt(SUM(P[ib,ie,im]^2 : ie\in 1..ne-1) + rdev[ib,im]^2)
#
  
  # Plot important data:
  plot_analysis_data \
    ( page, vset, masks, bands, illum_tb, wlen_tb, \
      savg, smin, smax, rdev, P
    )
  
def plot_analysis_data\
  ( page, vset, masks, bands, illum_tb, wlen_tb, \
    savg, smin, smax, rdev, P
  ):
  # Writes a file suitable for creating a plot of the results
  # {savg,smin,smax,rdev,P} from the analysis for page {page}, analysis run
  # {vset}, and a list of {nb} spectral bands with names
  # {bands[0..nb-1]}, comparing the different masks {masks[0..nm-1]} See
  # the top of this module for the format of this file.
  # 
  # Then calls the script {plot_joint_spectr.sh} to create a plot of that
  # spectrum as a ".png" file.
  #
  nb, ne, nm = numpy.shape(P)
  assert numpy.shape(bands) == (nb,), "bad {bands}"
  assert numpy.shape(savg) == (nb,nm,), "bad {savg}"
  assert numpy.shape(rdev) == (nb,nm,), "bad {rdev}"
  
  ddir = f"pages/{page}/vsets/{vset}"
  bash(f"mkdir -p {ddir}")
  dfile = f"{ddir}/plot.txt"
  with open(dfile, "w") as wr:
    wr.write(f"nb = {nb}\n")
    wr.write(f"nm = {nm}\n")

    wr.write("masks =")
    for im in range(nm): wr.write(f" {masks[im]}")
    wr.write("\n")

    for ib in range(nb):
      wr.write("%-15s" % bands[ib])
      illum = illum_tb[bands[ib]]
      wr.write(f" {illum:d}")
      wlen = wlen_tb[bands[ib]]
      wr.write(" %7.1f" % wlen)
      wr.write("  ")
      
      for im in range(nm):
        sum2_P = 0;
        for ke in range(1,ne): sum2_P += P[ib,ke,im]**2
        tdev = sqrt(sum2_P + rdev[ib,im]**2)

        wr.write("  ")
        wr.write(" %7.4f" % savg[ib,im])
        wr.write(" %7.4f" % smin[ib,im])
        wr.write(" %7.4f" % smax[ib,im])
        wr.write(" %7.4f" % tdev)
        wr.write(" %+7.4f" % P[ib,0,im])
      wr.write("\n")
    wr.close()
    
    bash(f"plot_analysis_data.sh {page} {vset}")
  return
  # ----------------------------------------------------------------------