-- $Id: point-01.lua,v 1.3 2006/06/06 19:41:43 dj Exp $
-- 
-- raygen specification for a single point source.  The point source
-- may have a single energy, multiple energies, a flat spectrum, or a
-- spectrum in a file.  This script does far too many things, but then
-- again, you won't have to write a separate script for all of them.
-- 
-- All options are passed to this script via the raygen override()
-- function which must set certain global variables which are
-- documented below.
-- 
-- There are two things that really must be specified: the output
-- spectrum and the source position.
-- 
-- Output Spectrum
-- ===============
-- 
--  This source can handle several different types of spectra.
--  This is all done by setting up a single variable, energy.
-- 
--    A single energy:
--    ----------------
--
--	Set energy equal to the energy [keV] and the variable `flux'
--	to the flux [photons/s/cm^2]
-- 
--        energy = 1.49
-- 
--   A list of discrete energies each with the same flux:
--   ----------------------------------------------------
--
--	Set energy equal to a table of energies [keV] and the variable
--	`flux' to the flux [photons/s/cm^2].
--
--	  energy = { 1.49, 4.51 }
-- 
--   A list of discrete energies with different fluxes:
--   --------------------------------------------------
--
--	Set energy equal to a table of tables, where each subtable 
--	contains the energy [keV] and flux [photons/s/cm^2]. e.g.
--
--	  energy = { {1.49, 0.5}, {4.51, 1.25} }
-- 
--   A list of discrete energies with different fluxes in a file
--   -----------------------------------------------------------
--
--	Set energy equal to a table, with a named element, picket,
--	which is itself a table with the following keys:
-- 
--	   - file	: name of RDB table
-- 	   - energy	: column which contains energy, in keV
--	   - flux	: column which contains flux
--	
--      the latter default to 'energy' and 'flux'.  if flux
--	is set to the string 'NONE', then all energies will have a
--      flux of 1
--
--	  energy = { picket = { file = 'foo.rdb' } }
-- 
--   A flat spectrum:
--   ----------------
-- 
--	Set energy to a table, with a named element, flat, which is
--	itself a table with the minimum and maximum energies [keV],
--	and the flux [photons/s/cm^2] (just like the raygen flat()
--	generator).  For example,
--
--	  energy = { flat = { 1.49, 5.0, 0.75 } }
-- 
--	The flux is optional and defaults to 1 photon/s/cm^2.
-- 
--   A spectrum in a file:
--   ---------------------
-- 
--	Set energy to a table, with a named element, specfile,
--	which is itself a table with the following keys:
-- 
--	   - file
--	   - format
-- 	   - emincol
--	   - emaxcol
--	   - fluxcol
--	   - units
--	   - scale
-- 
--	These have the same meaning as the arguments to the raygen
--	provided spectrum() function.  The file key must be provided.
--	For example,
-- 
--	  energy = { specfile = { file = 'foo.rdb' } }
-- 
-- 
--  You may combine all of the table based options:
-- 
--	energy = { 1.49 ; specfile = { ... }, flat = { 1.59, 5.0 } }
-- 
--  The different spectra are added together.

-- 
-- Source Position
-- ===============
-- 
--  The axial distance to the source is given by the variable `z',
--  where z is 0 at CAP datum A, and z increases from the source to
--  the focal plane.  If z is not specified, it is -Infinity.  Note
--  that the zero point is *not* at the telescope node.  This is an
--  issue only if the source is *really* close to the telescope.
-- 
--  The focal plane position of the source may be specified in one of
--  several coordinate systems.  The coordinate system is specified
--  via the `coord' variable, which may have one of the following
--  values:
--
--      raygen	- The native raygen polar coordinate system. This is the
--		  default.  Set the following variables:
-- 
--			theta - off axis angle in arcminutes
--			phi   - azimuthal angle in degrees
--
-- 
--	elaz	- OSAC elevation and azimuth. Set the following
--		  variables:
-- 
--			el - elevation in arcminutes
--			az - azimuth in arcminutes
-- 
--	msc	- Mirror Spherical Coordinates.  This is the on-orbit
--		  polar coordinate system. This makes an assumption
--		  about the relative orientations of the MSC and OSAC
--		  coordinate systems.  Set the following variables:
-- 
--			theta - off axis angle in arcminutes
--			phi   - azimuthal angle in degrees
--
--		   
--	offset	- orbit offset pointing. This also makes assumptions
--		  about the relative orientations of the orbit HRMA
--		  and OSAC coordinate systems. See the POG for more
--		  information on offset pointing. Specify the
--		  following variables:
--		  
--			yoff - Y offset in arcminutes
--			zoff - Z offset in arcminutes

$debug
AXAF_ROOT = getenv( 'AXAF_ROOT' )

dofile( AXAF_ROOT .. '/simul/lib/lua/RDB.lua' )
dofile( AXAF_ROOT .. '/simul/lib/lua/AXAFCoords.lua' )
dofile( AXAF_ROOT .. '/simul/lib/lua/store.lua' )


pfx = 'point-01.lua: '

flux = 1
coord = 'raygen'

override( )


sgen_idx = 0

function sgen_name( )
  sgen_idx = sgen_idx + 1
  return tostring( sgen_idx )
end

function gen_flat( table )

  if ( 'table' ~= type(table) ) then
    error( pfx .. "flat must be a table" )
  end 

  local emin = table[1]
  local emax = table[2]
  local _flux = table[3] or flux

  if ( nil == emin or 'number' ~= type(emin ) ) then
    error( pfx .. "flat spectrum has no or unparseable emin" )
  end

  if ( nil == emax or 'number' ~= type(emax ) ) then
    error( pfx .. "flat spectrum has no or unparseable emax" )
  end

  if ( 'number' ~= type(flux) ) then
    error( pfx .. "flat spectrum has unparseable flux" )
  end

  flat( sgen_name(), emin, emax, flux )

end

function picket( file, ecol, fluxcol )

  if ( nil == ecol ) then
     ecol = 'energy'
   end

  if ( nil == fluxcol ) then
     fluxcol = 'flux'
   end

   rdb, err = RDB:new( file )

   if ( nil == rdb ) then
      error( pfx .. err )
   end

   local cols

   if ( 'NONE' == fluxcol ) then
      cols = { ecol }
   else
      cols = { ecol, fluxcol }
   end

   local badcol = rdb:chk_cols( cols )
   if ( nil ~= badcol ) then
     error( file .. ': missing column ' .. badcol )
   end

   data = rdb:read()
   while data do
      if ( 'NONE' == fluxcol ) then
        flux = 1
      else
	flux = data[fluxcol]
      end
      mono( sgen_name(), data[ecol], flux )
      data = rdb:read()
   end

end


function energy_tbl_parse( idx, value )

  -- if the index is a number, then we've got a discrete energy
  if     ( 'number' == type(idx) ) then

    -- if the value is a number, it's the energy
    if     ( 'number' == type( value ) ) then
      mono( sgen_name(), value, flux )

    -- if the value is a table, should have { energy, flux } combo
    elseif ( 'table' == type( value ) ) then

      local _energy = value[1]
      local _flux = value[2] or flux

      if ( nil == _energy or 'number' ~= type(_energy ) ) then
         error( pfx .. "table element[" .. idx .. 
	           "] has no or unparseable energy" )
      end

      mono( sgen_name(), _energy, _flux )

    else
      error( pfx .. "unexpected value in energy array: " .. value )

    end

  -- if it's a string, must be one of the keywords we know
  elseif ( 'string' == type(idx) ) then

    -- flat spectrum
    if     ( 'flat' == idx ) then
      gen_flat( value )

    -- spectrum from file; just pass the table directly to spectrum()
    elseif ( 'specfile' == idx ) then
      spectrum( sgen_name(), value )

    elseif  ( 'picket' == idx ) then

      picket( value.file, value.energy, value.flux )    

    else
      error( pfx .. "unknown key in energy array: " .. idx )
    end

  end

end


function gen_spectrum()

--
  if ( "nil" == type( energy ) ) then

    error( pfx .. "energy not specified" )

--
  elseif ( "number" == type( energy ) ) then

    mono( 'mono', energy, flux )

  elseif ( "table" == type( energy ) ) then

    foreach( energy, energy_tbl_parse )

  end

end

function gen_pos()

  local attrs = {}

  if ( nil ~= z ) then
    attrs['z'] = z
  end

  coord = strlower(coord)

  if     ( 'raygen' == coord ) then

    attrs['theta'] = ( theta or 0 ) .. 'arcmin'
    attrs['phi']   = ( phi or 0 )  .. 'deg'

  elseif ( 'elaz' == coord ) then

    local el = ( el or 0 ) * min2rad
    local az = ( az or 0 ) * min2rad
    local theta, phi = osac_elaz_2_osac_polar( el, az )
    attrs['theta'], attrs['phi'] = osac_polar_2_raygen_polar( theta, phi )

  elseif ( 'msc' == coord ) then

    local theta = ( theta or 0 ) * min2rad
    local phi   = ( phi or 0 )   * deg2rad
    attrs['theta'], attrs['phi'] = MSC_to_raygen_polar( theta, phi )

  elseif ( 'offset' == coord ) then

    local el = - ( zoff or 0 ) * min2rad
    local az =   ( yoff or 0 ) * min2rad

    attrs['theta'], attrs['phi'] = raygen_elaz_2_raygen_polar( el, az )

  end


  point( "source", attrs )

end


function source( )

  begin_source( "point" )

  gen_pos()

  gen_spectrum( )

  end_source( "point" )
end