########################################
## spectline_subcont.g
##
## SUBTRACTING CONTINUUM FROM SPECTRAL LINE DATA CUBE 
##
## JHIBBARD 5/14/02
## called by specline_doall.g
######################################### 
##
## REQUIRES "BASENAME" AND "ROBUST" TO BE DEFINED FROM WITHOUT
## AND CUBE_MAPNAME HAS TO AGREE WITH WHAT IS IN specline_mapcube

## User Specified Parameters

AXIS			:= 4;              ## fequency axis 
                                           ## (get from summary on linecube)

ICHANSEL		:= [3,15,112,122]; ## continnum channels
                                           ## start1, stop1, start2, stop2

## Set up names for input and output images

CUBE_MAPNAME		:= spaste(BASENAME,'.R',ROBUST,'.line');
NOCONT_MAPNAME		:= spaste(BASENAME,'.R',ROBUST,'.nocont');
CONT_MAPNAME		:= spaste(BASENAME,'.R',ROBUST,'.cont');

######################################### 
# This script should be run several times, using different ICHANSEL
# to determin the channels that are really line-free. So we specify
# overwrite to be true 

# create image tool attached to line cube

include 'image.g';
im:=image(infile=CUBE_MAPNAME);

# Define region of image which is continuum, as defined by ICHANSEL.
#     Note that a union of regions is a compound region that
#     has to be made from two *world* regions, not pixel regions,
#     so we have to construct it from regions made with the
#     wrange function from the default regionmanager (drm).
#     See the reference manual on regionmanager for more on this.

# Store the image coordinate system so that regions 
# can be made in world coordinates

cs:=im.coordsys();                       

# Make first region cont1 for lower continuum channels

range1:=dq.quantity([ICHANSEL[1],ICHANSEL[2]],'pix');
cont1:=drm.wrange(pixelaxis=AXIS,range=range1,csys=cs);

# Make second region cont2 for upper continuum channels

range2:=dq.quantity([ICHANSEL[3],ICHANSEL[4]],'pix');
cont2:=drm.wrange(pixelaxis=AXIS,range=range2,csys=cs);


# Make region cont3 that specifies the union of the upper and 
# lower continuum channels

cont3:=drm.union(cont1,cont2);

# Form integrated continuum as moment=-1 of the line cube
# only using the continuum channels

im.moments(outfile=CONT_MAPNAME,
           moments=-1,
           axis=AXIS,
           region=cont3);

# Subtract continuum from line cube, writing continuum subtracted 
# cube to disk:

nocont:=imagecalc(outfile=NOCONT_MAPNAME, pixels=spaste('$im -',CONT_MAPNAME));

# View final image

nocont.view();

im.done();
cs.done();
nocont.done();

# Take a spectra through to check:

include 'imageprofilefitter.g';
pf := imageprofilefitter(infile=NOCONT_MAPNAME, infile2=CONT_MAPNAME);