#!/usr/bin/env julia

# Plot the 0th and 1st moment maps from the synthesized image.

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "--config"
    help = "a YAML configuration file"
    default = "config.yaml"
end

parsed_args = parse_args(ARGS, s)

import YAML
config = YAML.load(open(parsed_args["config"]))

using DiskJockey.constants
using DiskJockey.image
using DiskJockey.model
using HDF5

import PyPlot.plt
using LaTeXStrings
import Images

species = config["species"]
transition = config["transition"]
lam0 = lam0s[species*transition]
model = config["model"]
pars = convert_dict(config["parameters"], config["model"])

function plot_moment(mom::image.ZerothMoment, fname="zeroth_moment.png")

    # Image needs to be flipped along RA dimension
    ext = (mom.ra[end], mom.ra[1], mom.dec[1], mom.dec[end])

    fig, ax = plt[:subplots](nrows=1, figsize=(8,8))
    ax[:set_xlabel](L"$\Delta \alpha$ ('')")
    ax[:set_ylabel](L"$\Delta \delta$ ('')")

    # Flip the frame over the RA axis for Sky convention
    frame = flipdim(mom.data, 2)

    ax[:imshow](frame, extent=ext, interpolation="none", origin="lower", cmap=plt[:get_cmap]("PuBu"))

    plt[:savefig](fname)

end

# Assume that image has already been generated by synthesize_model.jl
im = imread()
skim = imToSky(im, pars.dpc)

# Do the velocity conversion here for plot labels
global nlam = length(skim.lams)
# convert wavelengths to velocities
global vels = c_kms * (skim.lams .- lam0)/lam0

moment = image.convert(image.ZerothMoment, skim)

plot_moment(moment)
