COLA recipes
What follows is a complete GIPSY task written in COLA. The first part is documentation. Further it contains information about how to read data from a text file, from the header of a set and from a GDS table. Also information can be found about the manipulation of strings. It writes GPLOT commands to a text file and after all the commands are written, GPLOT is started to read the commands from the file just created. At the end of the program, a system call to start an external PostScript viewer is described.
!
! COPYRIGHT (c) 1997
! Kapteyn Astronomical Institute
! University of Groningen, The Netherlands
! All Rights Reserved.
!
!
!#> atlas.dc1
!
!Program: atlas
!
!Purpose: Plot panels with data of a galaxy with given UGC number.
!
!Category: PLOTTING
!
!File: atlas.col
!
!Author: M.G.R. Vogelaar
!
!Keywords:
!
!
!**LOGFILE= Do you want messages in your log file? [Y]/N
!
! Default a lot of messages will be displayed in your
! log file. Most messages are generated by other
! applications that are started by this COLA file.
! You can suppress these messages if you type
! WHISPLOT LOGFILE=NO
!
!
! UGCNR= Enter UGC number: [Abort program]
!
! Enter the UGC number of the galaxy of which you want
! a panel plot. These UGC numbers are asked in a loop.
! The loop is aborted by pressing carriage return.
! Pending graphics of the first page half will not be
! written if NEWPAGE=N (default) was entered.
! See 'Description' to see which sets & files are needed
! to make a complete panel of plots.
!
!
! STANDARD= Do you want a standard UGC plot? [Y]/N
!
! Default is a half page mosaic with 6 panels containing
! the characteristics of a galaxy. The alternative
! (STANDARD=N) is a half page overview for the UGC number.
! Overview mosaics are usually plotted for an UGC number
! that is associated with more than one PGC number.
! the PGC number part in the default file names are
! replaced by 'OVERV'.
!
!
! PGCOK= Plot PGC .....? [Y]/N
!
! For each UGC number there is a text file on disk
! (uXXXX-galaxies.list) which contains one or more
! PGC numbers. A mosaic is made for each PGC number that
! enters the list with PGCOK=Y
!
! NEWPAGE= New page after first page half? Y/[N]
!
! Force an output with only one page half filled.
!
!
! GPLOTFILE= GPLOT input file: [Generated by program]
!
! This COLA program generates a so called 'input file'
! for GPLOT with all the necessary plot commands. After
! a name is entered, the program starts GPLOT which
! reads its commands from this file and writes its output
! to a PostScript file with the name entered in PSFILE=
!
!
! OVERWRITE= File [GP ..... .txt] exists. Overwrite? Y/[N]
!
! This prompt is generated if you specified a file name
! in GPLOTFILE= that already exists.
!
!
! PSFILE= PostScript output file: [Generated by program]
!
!
! GRDEVICE= Enter GPLOT device: [PCPSFILE]
!
! The destination of the plot. Input is the same as
! other GRDEVICE= keywords in GIPSY. If you write
! to a PostScript file( e.g. PCPSFILE), you can
! view the result (using Ghostview or Ghostscript) before
! sending it to a printer. Note that the optical image is
! a colour image, so the device must be a colour device.
! However a colour PostScript file can be send to a
! black and white printer.
!
!
!
!
!Description: The program plots a mosaic of plots that contain information
! about a galaxy. Per galaxy there are two rows of three plots
! and there are panels of two galaxies on one page.
! There is a loop over UGC numbers (UGCNR=). For each UGC number
! there is one, or more than one, PCG (Paturel Galaxy Catalog)
! number. For each PCG number a mosaic will be plotted.
! The PCG numbers are extracted from an external ASCII file.
! The name of that file is always "uXXXX-galaxies.list"
! where XXXX is the UGC number. For each PGC number the
! user is asked to confirm or reject the entry PGCOK=.
! If the list with PGC numbers does not exist, the program
! is aborted.
! Two PCG numbers will fill a page. For each (PostScript) page
! a name is created for the PostScript file (name starts
! with PS) and for the GPLOT input file (starts with GP).
! The rest of the name contains the UGC and PGC numbers of
! the objects that are plotted. An example of a GPLOT file
! name is 'GPu9211p51357u9210p51358.txt'. The defaults
! can be changed with keywords GPLOTFILE= and PSFILE=.
! Usually the mosaics
! contain standard plots like the optical image , a HI map,
! a velocity field etc., but the user can also select an
! overview of an UGC number (STANDARD=N). This is also a
! mosaic which occupies half a page, but another set of
! keywords is asked to create the panels.
!
! The standard page half consists of 6 plots. All plots
! need a source name. Sometimes this is the name of a
! GIPSY set, sometimes it is an ASCII file. First the
! program checks which files are available. In the next table
! a listing is made which files are necessary:
!
! First an UGC
! number is required. With this number all the necessary
! filenames are generated. Below is a list with the plot
! contents ('XXXX' is the UGC number and YYYYY is the PGC
! number). Only the HI image is required. The other data
! is optional.
!
! row col filename/set description
! ================================================================
! 1 1 uXXXXopt optical image
! or: .optical/ugcXXXX.fits
! 1 2 uXXXXs30-2dim 0 Total HI, 30 arcsec, col. map
! and: uXXXXs30-2dim 8 SNR levels for mask
! 1 3 uXXXXs30-2dim 1 I.W.M vel, 30 arc., col
! 2 1 fluxXXXXs60han.dat global profile
! 2 2 uXXXXs30-2dim 9 continuum
! 2 3 uXXXXs30xvYYYYY 0 XV map
!
! The box for the optical/HI/I.W.M vel and continuum images
! is the same. This box is set only once and stored in the
! header of the HI set.
!
! The colour settings can be adjusted separately for both
! total HI data and the i.w.m velocity.
!
! The optical image is obtained from a Digitized Sky Survey
! FITS file. This file must be in a subdirectory 'optical'
! and has name ugcXXXX.fits (e.g.ugc9211.fits).
!
! This is an example of the contents of a so called
! 'info file', e.g. u9211p51358info.txt.
! The file is created by some other program and contains
! the title of the plot (displayed at the top left), the
! beam sizes, the channel separation (used to plot a beam
! in the XV plot).
!
!
! !-----------------------INFO FILE----------------------------
! !
! ! Example of an info file. Comments must start with a '!' as
! ! the FIRST character of a line
! !
! !------------------------------------------------------------
! ! Name for title at top of plot
! TITLE=UGC 9211
! ! Beam minor axis in ARCSEC
! BMMINFULL=11.5
! ! Major axis
! BMMAJFULL=16.1
! ! For 30'' resolution
! BMMIN30=29.7
! BMMAJ30=30.1
! ! Channel separation in KM/S
! CHSEP=4.15
! ! Ellipse parameters: inclination, position angle and major axis.
! ! Major axis of overlay ellipse always in ARCMIN
! ELLMAJ=1.0
! ELLPA=30.0
! ! Inclination determines length of minor axis
! ELLINC=40.0
! ! Optical centre used to centre ellipse in continuum plot
! ! Note that this position follows the GIPSY convention for positions
! RAOPT=* 15 41 22
! DECOPT=* 67 24 30
! GALTYPE=SC
! ! System velocity always in KM/S
! VSYS=425.0
!
!
! The colours of the HI map at 30'' is set by GIDS unless
! colours were stored in the header of this set (at set level).
! The HI 30'' set is modified by transferring all data
! greater then a certain signal to noise level. The snr levels
! are stored in subset 8 of the 2-dim map.
!
! The global profile data is read from a file named:
! fluxuXXXXs60hanpYYYYY.dat. If the PGC number is not available,
! the name without the p-extension is used. The data in this
! file is in the format x, y. The data in the y array is
! multiplied with 5.0 to scale from W.U. to mJy's.
!
! The XV map is named uXXXXs30xvYYYYY. Note that the character
! 'p' is omitted. The map must be extracted at DEC equal to 0.
!
!Notes:
!
!Example:
!
!Updates: Jun 6, 1997: VOG, Document created
!
!#<
string setname(9) !
string tempname
string plotid(9) ! Identification for use in prompt
string ugcnr ! UGC number as strings
string pgcnr(20) ! Max. number of PGC nr's per UGC is 20
string gplotfile ! Name of GPLOT Ascii file
string galaxylist ! Name of list with galaxy info
string infofile ! Plot text from this file below plots
string line ! One line from galaxylist file
string progname ! Name of this COLA file
string set ! setname without subset info
string boxstr ! Box info from table in header
string plottitle(9) ! Title above each plot
string toptitle ! Title above panel
string grdevice ! Name of GPLOT plot device
string axdeltaX ! Value and units of delta in major ticks
string axdeltaY
string axposX ! Position of first label along X
string axposY
string axformX ! Format for axis labels
string axformY
string axtitleX ! Label title
string axtitleY
string*20 offunits ! Units along XV's offset axis
string*20 velunits
string*20 unitsdelta
string*20 unitspos
string defstr ! A string for defaults
string tmpset ! Temp. set to store SNR correlated IWMV/tHI data
string basename ! Base name for temp. sets
string clipstr ! Default clip for HI plots
string psfile ! Name of PostScript result
string infoitem(10) ! Names of info items from info file
string infocont(10) ! Text of those items
string substr ! Dummy part of a string
string keyword ! Keyword from info file
string fitsfile ! get optical image from fits file
string psviewer
string*256 copyline ! Copy (long) lines from temp file to GPLOT file
string levels
string galtype
string xrange, yrange
string*20 infokey(10)
string yp
logical exist(9) ! Can set/files be opened?
logical gplotfileexist
logical galaxylistexist ! Do these files exist?
logical fromheader ! Is a (header) table available?
logical foundPGC ! Is a PGC number found in Ascii file?
logical foundkey ! found a matching keyword in info file
logical morelines ! Are there more lines to read?
logical lg, fileexist, ok, overwrite ! result values
logical pgcok ! Want to plot this PCG?
logical viewps ! Start ghostview?
logical logfile ! Display log file messages?
logical newpage ! New page after one galaxy?
logical stopprogram
logical standard
logical plotmade
logical first
integer tempfile ! temporary text file for GPLOT input
integer datfile
integer r ! result value
integer gpfp ! Pointer to GPLOT Ascii file
integer galfp ! Pointer to galaxy list file
integer infofp ! Pointer to info file
integer subsets(1024) ! Max. number of input subsets
integer nsubs ! Number of input subsets
integer setlevel ! Top header level
integer pn ! One dim. plot number
integer pp ! 1-dim plot position (as pn)
integer subnr ! A subset
integer gplotoutput ! The log file info destination
integer i ! counter
integer numlines ! info lines read from 'info' file
integer len1, len2 ! Length of strings
integer maxinfolines ! Max. number of lines from info file
integer numPGC
integer pagehalf ! Two panel plots on one page
integer panelnr
integer col, row
integer ch ! Character counter
integer n
integer numkeys
integer keynr
real x, y ! positions in plot n mm
real xoff
real xstart, ystart ! Start pos. of first plot in mm
real xgap, ygap ! Space between plots in mm
real Xsizeplot, Ysizeplot ! Width and height of panels in mm
real chtoplabel, chlabels ! Character heights in mm
real chplottitle ! Titles above plot
real chinfo ! Char. height of info text
real chinfospace ! Space in mm between these lines
real charrow
real xmargin, ymargin ! Space between plot and frame
real titlegap ! Space between frame and plot title
real colbarheight ! Height of colour bar in mm
real colbarwidth ! Width of colour bar in mm
real snr
real snrHI ! Threshold for total HI maps
real snrIWM ! Threshold for IWM vel maps
real datamin ! Default values for clip in HI plots
real datamax
real fctor ! scaling factor for global profile
real beamDx, beamDy ! Number between 0 and 1 which sets
! the offset of the center of a beam
! in terms of 'xsize' and 'ysize'
! wrt. the lower left corner of a plot
real A4Xsize ! Sizes of A4 in mm
real A4Ysize
real A4Xmargin ! Bottom and top margin for a page half
real A4Ymargin
real plotdeltaX ! Variable used yo center plots
real plotdeltaY
real levstep
real numlevels
real xdatmax , xdatmin
real ydatmax , ydatmin
real xdat , ydat
double bmmin ! Minor axis beam
double bmmaj ! Major
double chsep ! Channel separation in km/s
double ellinc
double ellpa
double ellmaj, ellmin
double vsys
double scale ! scaling factor from wu to nHI
double cdelt1, cdelt2
double crval1, crval2
double raopt, decopt
double keycon(9)
!--------------------------------------------------
! Miscellaneous constants
!--------------------------------------------------
progname = "ATLAS" ! Name of this program
setlevel = 0 ! For header items on top level
snr = 3.0
snrIWM = 3.0 ! Default clip for Intensity weighted mean velocity
snrHI = 3.0 ! Default clip for Total HI maps
tempname = "dummygplot.txt"
A4Ysize = 296.5
A4Xsize = 209.5
A4Xmargin = 25.0
A4Ymargin = 8.0
xmargin = 1.0
ymargin = 1.0
Xsizeplot = 45.0
Ysizeplot = 45.0
plotdeltaX = 15.0
plotdeltaY = 12.0
titlegap = 2.0
chtoplabel = 4.0 ! Character heights in mm
chlabels = 2.0 ! Physical coordinate labels along the axes
chplottitle = 2.2 ! Titles above each plot
chinfo = 3.3 ! Text character height info file lines
chinfospace = 0.9 ! Space in mm between lines from info text file
charrow = 8.0
plotid(1) = "Opt."
plotid(2) = "HI30''"
plotid(3) = "Vel.field"
plotid(4) = "Glob.prof."
plotid(5) = "Continuum"
plotid(6) = "XV"
plottitle(1) = "Optical"
plottitle(2) = "Total HI (30 arcsec res.)"
plottitle(3) = "I.w.m. velocity (30 arcsec res.)"
plottitle(4) = "Global HI profile"
plottitle(5) = "Continuum (30 arcsec res.)"
plottitle(6) = "Pos.-vel. plot along major axis"
toptitle = ""
psviewer = "gv"
beamDx = 0.9
beamDy = 0.1
numlevels = 10
infokey(1) = "BMMIN30="
infokey(2) = "BMMAJ30="
infokey(3) = "RAOPT="
infokey(4) = "DECOPT="
infokey(5) = "ELLINC= "
infokey(6) = "ELLPA="
infokey(7) = "ELLMAJ="
infokey(8) = "VSYS="
infokey(9) = "CHSEP="
numkeys = 9
!--------------------------------------------------
! Set output level for log file info
!--------------------------------------------------
r = setdef( 2 ) ! Set default level to hidden
logfile = false
read "Do you want GPLOT output messages? Y/[N]" logfile
if (logfile) then
gplotoutput = 3
else
gplotoutput = 0
cif
r = setdef( 1 )
pagehalf = 1
repeat
!--------------------------------------------------
! File names that are used to plot data all have
! an UGC (or PGC) number in it. Start the program by
! asking this UGC number.
!--------------------------------------------------
ugcnr = ""
read "Enter UGC number: [Abort program]" ugcnr
cancel ugcnr
if (len(ugcnr) = 0) then
halt
cif
!-------------------------
! Does HI-30'' set exists?
!-------------------------
pn = 2
setname(pn) = "u%ugcnr_s30-2dim"
exist(pn) = exists(setname(pn))
if (not exist(pn)) then
write "[%progname] Set %setname(pn) cannot be opened! Aborting program..."
halt
else
pn = 3
setname(pn) = setname(2)
exist(pn) = exists(setname(pn))
pn = 5
setname(pn) = setname(2)
exist(pn) = exists(setname(pn))
cif
!-------------------------
! Does optical set exists?
!-------------------------
pn = 1
setname(pn) = "u%ugcnr_opt" ! A 2-dim set
exist(pn) = exists( setname(pn) )
if (not exist(pn)) then
write "[%progname] Could not find set [%setname(pn)], start program DSSLOAD"
fitsfile = "./optical/ugc%ugcnr.fits"
if (not existf(fitsfile)) then
write "[%progname] Could not find fits file [%fitsfile]"
else
tmpset = "u%ugcnr_opttmp"
"DSSLOAD FILENAME=%fitsfile OUTSET=%tmpset"
! reproject to sky system of HI-30''
r = rhead( "%setname(2)", 0, "CRVAL1", crval1 )
r = rhead( "%setname(2)", 0, "CRVAL2", crval2 )
"REPROJ INSET=%tmpset BOX= DEFSET= ROTATEONLY=N
OUTPOS=U %crval1 U %crval2 CDELT=
SKYSYS=1 EPOCH=1950 PROSYS=8 CROTA=0.0 OUTBOX=
SPEEDMAT= DATAMODE= OUTSET=u%ugcnr_opt"
"DELETE INSET=%tmpset; OK=Y;"
exist(pn) = exists( setname(pn) )
if (exist(pn)) then
nsubs = getsas( "%setname(pn)", set, subsets )
exist(pn) = (nsubs > 0)
if (not exist(pn)) then
write "[%progname] Something wrong with set [%setname(pn)]"
cif
cif
cif
cif
standard = true
read "Do you want a standard plot? [Y]/N" standard
cancel standard
!--------------------------------------------------
! Open file uXXXX-galaxies.list (XXXX=is the actual
! UGC number) to extract PGC number(s).
! For each PGC entry that is found in the file,
! the user can reject the number. The number of
! PGC numbers determine how many (sub) loops
! there are for one UGC number. Check the existence
! of the 'galaxy list'.
!--------------------------------------------------
galaxylist = "u%ugcnr_-galaxies.list"
galaxylistexist = existf( galaxylist )
if (not galaxylistexist) then
write "[%progname] Could not open galaxy list [%galaxylist]"
write "therefore no PGC number(s) could be found. Aborting program!"
halt
cif
if (not standard) then
!--------------------------------------------------
! There are no PGC numbers from a galaxy list but
! an overview is wanted. Substitute for this
! overview the PGC number 'OVERV' and the number
! of PGC numbers 1.
!--------------------------------------------------
numPGC = 1
pgcnr(1) = "OVERV"
else
!--------------------------------------------------
! The file is found, now read the lines and extract
! the PGC number(s).
!--------------------------------------------------
i = 1
galfp = fopen( galaxylist, "r" )
repeat
pgcnr(i) = ""
r = fread( galfp, line )
morelines = (r >= 0)
if (morelines) then
foundPGC = ( "PGC" = substr1(line, 1) )
if (foundPGC) then
pgcnr(i) = substr1(line, 2)
pgcok = true
read "Plot PGC %pgcnr(i)? [Y]/N" pgcok
cancel pgcok
if (pgcok) then
i = i + 1
else
pgcnr(i) = ""
cif
cif
cif
until (not morelines)
r = fclose( galfp ) ! Close galaxy list
numPGC = i - 1 ! The number of selected PGC galaxies
!--------------------------------------------------
! Generate some output for control
!--------------------------------------------------
if (numPGC < 1) then
write "[%progname] Could not find any PGC numbers. The program is aborted!"
halt
else
write "[%progname] You selected the following PGC number(s):"
for i = 1, numPGC
write "%pgcnr(i)"
cfor
cif
cif ! End PGC loop == not an overview
!--------------------------------------------------
! Start a loop over all the PGC numbers. For each
! number a number of keywords must be specified
! before a 'half page' plot is made.
!--------------------------------------------------
for i = 1, numPGC
if (pagehalf = 1) then
!--------------------------------------------------
! Open a temporary file. Write all GPLOT commands
! to this file. Rename it later with unique name.
!--------------------------------------------------
tempfile = fopen( tempname, "w" )
if (tempfile < 0) then ! Can file be opened?
write "[%progname] Cannot open dummy text file! Aborting program..."
halt
cif
cif
if (not standard) then
!--------------------------------------------------
! Plot a special overview panel. Note that for an
! overview: numPGC=1
!--------------------------------------------------
else
!--------------------------------------------------
! Before we can plot any frames, we have to know
! what the box is. The box is the same for all the
! plots with spatial axes so we have to set it
! only once. This is done by displaying the total
! HI map (at 30'') and starting program
! MINBOX -the map is usually the result of MOMENTS
! and therefore the subset number for this map is 0.
! Then the box is stored as a string
! in a table in the descriptor of this HI set.
! If this COLA file encounters a box description
! in its header, then it will prompt with this value
! as default.
!--------------------------------------------------
subnr = 0
boxstr = "-10 -10 10 10"
pn = 2
nsubs = getsas( "%setname(pn) %subnr", set, subsets )
if (not (nsubs > 0)) then
write "[%progname] Something wrong with set [%setname(pn) %subnr]"
else
fromheader = (existt( setname(pn), subsets(1), "MBOX") > 0)
boxstr = ""
if (fromheader) then
r = rcolumn( setname(pn), subsets(1), "MBOX", "MINBOX", boxstr, 1, 1 )
fromheader = (len(boxstr) > 0)
if (fromheader) then
read "%Plotid(2): Retrieve BOX from this header? [Y]/N" fromheader
cancel fromheader
cif
cif
if (not fromheader) then
!--------------------------------------------------
! Let MINBOX write a box string in this header.
!--------------------------------------------------
"VIEW INSET=%setname(pn) %subnr"
"MINBOX INSET=%setname(pn) %subnr OUTSET= SQUARE=yes"
r = rcolumn( setname(pn), subsets(1), "MBOX", "MINBOX", boxstr, 1, 1 )
cif
! Purify this string
boxstr = substr2( boxstr, 1, len(boxstr) )
!--------------------------------------------------
! String 'boxstr' contains description of current
! box, and is a global variable for a couple of
! panel plots.
!--------------------------------------------------
cif
!--------------------------------------------------
! We have defined a box now and want to plot frames
! around the plots where both axes are spatial.
! There is a default labeling, but usually the
! user wants to adjust this. If a frame is set up,
! the characteristics are stored in the header of
! the 30 '' HI set.
!--------------------------------------------------
axdeltaX = ""
axdeltaY = ""
axposX = ""
axposY = ""
axformX = ""
axformY = ""
axtitleX = ""
axtitleY = ""
if (exist(pn)) then
fromheader = (existt(setname(pn), subsets(1), "AXES") > 0)
if (fromheader) then
ok = true
read "%Plotid(pn): Retrieve AXIS properties from header?: [Y]/N" ok
cancel ok
if (ok) then
r = rcolumn( setname(pn), subsets(1), "AXES", "AXDELTX", axdeltaX, 1, 1 )
r = rcolumn( setname(pn), subsets(1), "AXES", "AXDELTY", axdeltaY, 1, 1 )
r = rcolumn( setname(pn), subsets(1), "AXES", "AXPOSX", axposX, 1, 1 )
r = rcolumn( setname(pn), subsets(1), "AXES", "AXPOSY", axposY, 1, 1 )
r = rcolumn( setname(pn), subsets(1), "AXES", "AXFORMX", axformX, 1, 1 )
r = rcolumn( setname(pn), subsets(1), "AXES", "AXFORMY", axformY, 1, 1 )
r = rcolumn( setname(pn), subsets(1), "AXES", "AXTITX", axtitleX, 1, 1 )
r = rcolumn( setname(pn), subsets(1), "AXES", "AXTITY", axtitleY, 1, 1 )
cif
cif
repeat
r = list( 0 )
"GPLOT COMMAND=device x11;inset %setname(pn) %subnr;box %boxstr;autoscale;
axdelta %axdeltaX;axpos %axposX;axform %axformX;axis bp;axtitle %axtitleX;
axdelta %axdeltaY;axpos %axposY;axform %axformY;axis lp;axtitle %axtitleY;axis t;axis r;quit"
r = list( 1 )
ok = true
read "Labels etc. ok? [Y]/N" ok
cancel ok
if (not ok) then
defstr = substr2(axdeltaX, 1, len(axdeltaX) )
read "Give delta(+units) for X-axis labeling: [%defstr]" axdeltaX
cancel axdeltaX
defstr = substr2(axposX, 1, len(axposX) )
read "Give position of 1th X-axis label: [%defstr]" axposX
cancel axposX
defstr = substr2(axformX, 1, len(axformX) )
read "Give format of X-axis labels: [%defstr]" axformX
cancel axformX
defstr = substr2(axtitleX, 1, len(axtitleX) )
read "Give title along X-axis: [%defstr]" axtitleX
cancel axtitleX
defstr = substr2(axdeltaY, 1, len(axdeltaY) )
read "Give delta(+units) for Y-axis labeling: [%defstr]" axdeltaY
cancel axdeltaY
defstr = substr2(axposY, 1, len(axposY) )
read "Give position of 1th Y-axis label: [%defstr]" axposY
cancel axposY
defstr = substr2(axformY, 1, len(axformY) )
read "Give format of Y-axis labels: [%defstr]" axformY
cancel axformY
defstr = substr2(axtitleY, 1, len(axtitleY) )
read "Give title along Y-axis: [%defstr]" axtitleY
cancel axtitleY
cif
until (ok)
!--------------------------------------------------
! Store these values in the header
!--------------------------------------------------
if (fromheader) then
r = deletet( setname(pn), subsets(1), "AXES" )
cif
r = len( axdeltaX )
r = createc( setname(pn), subsets(1), "AXES", "AXDELTX", "CHAR%r",
"Space between major tickmarks in X", " " )
r = wcolumn( setname(pn), subsets(1), "AXES", "AXDELTX", axdeltaX, 1, 1 )
r = len( axposX )
r = createc( setname(pn), subsets(1), "AXES", "AXPOSX", "CHAR%r",
"Position of first (complete) label in X", " " )
r = wcolumn( setname(pn), subsets(1), "AXES", "AXPOSX", axposX, 1, 1 )
r = len( axformX )
r = createc( setname(pn), subsets(1), "AXES", "AXFORMX", "CHAR%r",
"Format of the labels in X", " " )
r = wcolumn( setname(pn), subsets(1), "AXES", "AXFORMX", axformX, 1, 1 )
r = len( axtitleX )
r = createc( setname(pn), subsets(1), "AXES", "AXTITX", "CHAR%r",
"Title along X-axis", " " )
r = wcolumn( setname(pn), subsets(1), "AXES", "AXTITX", axtitleX, 1, 1 )
r = len( axdeltaY )
r = createc( setname(pn), subsets(1), "AXES", "AXDELTY", "CHAR%r",
"Space between major tickmarks in Y", " " )
r = wcolumn( setname(pn), subsets(1), "AXES", "AXDELTY", axdeltaY, 1, 1 )
r = len( axposY )
r = createc( setname(pn), subsets(1), "AXES", "AXPOSY", "CHAR%r",
"Position of first (complete) label in Y", " " )
r = wcolumn( setname(pn), subsets(1), "AXES", "AXPOSY", axposY, 1, 1 )
r = len( axformY )
r = createc( setname(pn), subsets(1), "AXES", "AXFORMY", "CHAR%r",
"Format of the labels in Y", " " )
r = wcolumn( setname(pn), subsets(1), "AXES", "AXFORMY", axformY, 1, 1 )
r = len( axtitleY )
r = createc( setname(pn), subsets(1), "AXES", "AXTITY", "CHAR%r",
"Title along Y-axis", " " )
r = wcolumn( setname(pn), subsets(1), "AXES", "AXTITY", axtitleY, 1, 1 )
cif
!--------------------------------------------------
! Open the info file and store its contents
! The info file (e.g. u9211p51358info.txt) contains
! both the UGC and the PGC number
!--------------------------------------------------
infofile = "u%ugcnr_p%pgcnr(i)_info.txt"
numlines = 0
toptitle= "?"
if ( NOT existf(infofile) ) then
write "[%progname] File [%infofile] does not exist. Defaults substituted!"
else
infofp = fopen( infofile, "r" )
write "[%progname] Reading additional information from file [%infofile]"
i = 1
repeat
r = fread( infofp, line )
morelines = (r >= 0)
if (morelines) then
foundkey = false
!--------------------------------------------------
! If a string starts with a '!' then it is comment.
!--------------------------------------------------
if (substr2(line, 1, 1) = "!") then
foundkey = true
else
len2 = len( line )
cif
if (not foundkey) then
keyword = "TITLE="
len1 = len( keyword )
substr = substr2( line, 1, len1 )
if (substr = keyword ) then
foundkey = true
toptitle = substr2( line, len1+1, len2-len1 )
cif
cif
if (not foundkey) then
keyword = "GALTYPE="
len1 = len( keyword )
substr = substr2( line, 1, len1 )
if (substr = keyword ) then
foundkey = true
galtype = substr2( line, len1+1, len2-len1 )
galtype = substr2( galtype, 1, len(galtype) )
cif
cif
! Read the numbers
if (not foundkey) then
keynr = 1
repeat
len1 = len( infokey(keynr) )
substr = substr2( line, 1, len1 )
if (substr = infokey(keynr) ) then
foundkey = true
substr = substr2(line, len1+1, len2-len1)
keycon(keynr) = atof( substr2(line, len1+1, len2-len1) )
cif
if (foundkey) then
write "%infokey(keynr)_%keycon(keynr)"
cif
keynr = keynr + 1
until (foundkey or (keynr > numkeys))
cif
if (foundkey) then
cif
cif
until (not morelines)
r = fclose( infofp )
numlines = i - 1
bmmin = keycon(1)
bmmaj = keycon(2)
raopt = keycon(3)
decopt = keycon(4)
ellinc = keycon(5)
ellpa = keycon(6)
ellmaj = keycon(7)
ellmin = ellmaj*COS(RAD(ellinc))
vsys = keycon(8)
chsep = keycon(9)
cif
!--------------------------------------------------
! Plot a label at the top: Overview galaxy .....
!--------------------------------------------------
! x = A4Xmargin + plotdeltaX+Xsizeplot + Xsizeplot/2.0
! y = A4Ymargin + plotdeltaY+2.0*Ysizeplot + plotdeltaY/2.0
! if (pagehalf = 1) then
! y = y + 2.0 * (plotdeltaY+Ysizeplot) + plotdeltaY/2.0
! cif
! read "Enter title of plot: [%toptitle_]" toptitle
! cancel toptitle
! r = fwrite( tempfile, "mm" )
! r = fwrite( tempfile, "move %x, %y" )
! r = fwrite( tempfile, "justifi 0.5" )
! r = fwrite( tempfile, "lwidth 2" )
! r = fwrite( tempfile, "charheight %chtoplabel" )
! r = fwrite( tempfile, "text %toptitle" )
! r = fwrite( tempfile, "justifi 0.0" )
! r = fwrite( tempfile, "lwidth 1" )
! r = fwrite( tempfile, "charheight %chlabels" )
! r = fwrite( tempfile, "world" )
!--------------------------------------------------
! Standard 'atlas' plot
!--------------------------------------------------
for panelnr = 1, 6
plotmade = false
if (panelnr < 4) then
row = 1
col = panelnr
else
row = 2
col = panelnr - 3
cif
x = A4Xmargin + (col-1.0) * (plotdeltaX+Xsizeplot)
y = A4Ymargin + (2.0-row) * (plotdeltaY+Ysizeplot)
if (pagehalf = 1) then
y = y + 2.0 * (plotdeltaY+Ysizeplot) + plotdeltaY/2.0
cif
r = fwrite( tempfile, "location %x %y" )
r = fwrite( tempfile, "charheight %chlabels" )
pn = panelnr
levels = ""
if (pn = 1 and exist(1) ) then
write "[%progname] Plot %pn is set [%setname(pn)], title: %plottitle(pn)"
r = fwrite( tempfile, "inset %setname(pn)" )
! Box for this optical image is the same (in phys. coords. as HI maps
r = fwrite( tempfile, "box %boxstr" )
r = fwrite( tempfile, "xsize %Xsizeplot" )
r = fwrite( tempfile, "ysize %Ysizeplot" )
fromheader = (existt(setname(pn), 0, "LUTS") > 3)
if (fromheader) then
read "%Plotid(pn): Retrieve colour lut from header? [Y]/N" fromheader
cancel fromheader
cif
if (not fromheader) then
"VIEW INSET=%setname(pn) BOX=%boxstr"
"GPLOT COMMAND=device null;inset %setname(pn);pause;getlut;storelut;quit"
cif
r = fwrite( tempfile, "getlut header" )
r = fwrite( tempfile, "colplot" )
r = fwrite( tempfile, "axdelta %axdeltaX" )
r = fwrite( tempfile, "axpos %axposX" )
r = fwrite( tempfile, "axformat %axformX" )
r = fwrite( tempfile, "axis bp" )
r = fwrite( tempfile, "axtitle %axtitleX" )
r = fwrite( tempfile, "axdelta %axdeltaY" )
r = fwrite( tempfile, "axpos %axposY" )
r = fwrite( tempfile, "axformat %axformY" )
r = fwrite( tempfile, "axis lp" )
r = fwrite( tempfile, "axtitle %axtitleY" )
r = fwrite( tempfile, "axis t" )
r = fwrite( tempfile, "axis r" )
r = fwrite( tempfile, "mm" )
r = fwrite( tempfile, "justifi 1.0" )
r = fwrite( tempfile, "move %x+%Ysizeplot %y+%Ysizeplot+%titlegap" )
r = fwrite( tempfile, "charheight %chplottitle" )
r = fwrite( tempfile, "text %plottitle(pn)" )
r = fwrite( tempfile, "justifi 0.0" )
r = fwrite( tempfile, "move %x %y+%Ysizeplot+%titlegap" )
r = fwrite( tempfile, "charheight %chtoplabel" )
r = fwrite( tempfile, "lwidth 2" )
r = fwrite( tempfile, "text %toptitle" )
r = fwrite( tempfile, "lwidth 1" )
r = fwrite( tempfile, "charheight %chlabels" )
r = fwrite( tempfile, "colour background")
r = fwrite( tempfile, "fastyle 1" ) ! solid
r = fwrite( tempfile, "move %x+2.0 %y+%Ysizeplot+%titlegap-5.0" )
r = fwrite( tempfile, "rectangle 15 3.5" )
r = fwrite( tempfile, "move %x+3.0 %y+%Ysizeplot+%titlegap-4.0" )
r = fwrite( tempfile, "colour foreground")
r = fwrite( tempfile, "text type: %galtype" )
r = fwrite( tempfile, "world" )
plotmade = true
cif ! panelnr one
if ( (pn = 2 or pn = 3 or pn = 5) and exist(2) ) then
!--------------------------------------------------
! The colours in the HI 30'' are set by GIDS. The
! colour properties are stored in the header of
! the total HI at 30'' resolution at subset level
! The I.W.M vel map is treated the same way.
!--------------------------------------------------
tmpset = "TMP%pn_u%ugcnr_p%pgcnr(i)"
if (pn = 2) then
subnr = 0
elsif (pn = 3) then
subnr = 1
else
subnr = 9
cif
write "[%progname] Plot %pn is set [%setname(pn)], title: %plottitle(pn)"
nsubs = getsas( "%setname(pn) %subnr", set, subsets )
read "%Plotid(pn): Enter SNR threshold for map: [%snr]" snr
cancel snr
r = list( 0 )
"MNMX INSET=%setname(pn) %subnr; BOX=%boxstr CHANGE=YES"
r = list( 1 )
r = rhead( set, subsets(1), "DATAMAX", datamax )
if (r < 0) then
datamax = 0.0
cif
r = rhead( set, subsets(1), "DATAMIN", datamin )
if (r < 0) then
datamin = 0.0
cif
r = rcolumn( setname(pn), subsets(1), "LEVTAB", "LEVELS", levels, 1, 1 )
fromheader = (len(levels) > 0)
if (not fromheader) then
levstep = (datamax - datamin) / numlevels
if (pn = 3) then
levels = "%.2datamin_:%.2datamax_:10"
else
levels = "%.2datamin_:%.2datamax_:%.2levstep"
cif
cif
levels = substr2( levels, 1, len(levels) )
read "%Plotid(pn): Enter levels for this plot [%levels]" levels
cancel levels
if (fromheader) then
r = deletet( setname(pn), subsets(1), "LEVTAB" )
cif
r = len( levels )
r = createc( setname(pn), subsets(1), "LEVTAB", "LEVELS", "CHAR%r",
"Levels in mapunits for contours", " " )
r = wcolumn( setname(pn), subsets(1), "LEVTAB", "LEVELS", levels, 1, 1 )
if (pn = 5) then
r = fwrite( tempfile, "task combin RESULT01=$1 RESULT02= SET01=%setname(pn) %subnr SET02=%setname(pn) 8 SETOUT01=%tmpset 1 BOX01=%boxstr BOX02=" )
else
r = fwrite( tempfile, "task combin RESULT01=ifgt($2,%snr,$1,blank) RESULT02= SET01=%setname(pn) %subnr SET02=%setname(pn) 8 SETOUT01=%tmpset 1 BOX01=%boxstr BOX02=" )
cif
axdeltaX = substr2( axdeltaX, 1, len(axdeltaX) )
axdeltaY = substr2( axdeltaY, 1, len(axdeltaY) )
axposX = substr2( axposX, 1, len(axposX) )
axposY = substr2( axposY, 1, len(axposY) )
axformX = substr2( axformX, 1, len(axformX) )
axformY = substr2( axformY, 1, len(axformY) )
axtitleX = substr2( axtitleX, 1, len(axtitleX) )
axtitleY = substr2( axtitleY, 1, len(axtitleY) )
r = fwrite( tempfile, "inset %tmpset 1" )
r = fwrite( tempfile, "box %boxstr" )
r = fwrite( tempfile, "xsize %Xsizeplot" )
r = fwrite( tempfile, "ysize %Ysizeplot" )
r = fwrite( tempfile, "levels %levels" )
r = fwrite( tempfile, "grayscale" )
r = fwrite( tempfile, "contours" )
! r = fwrite( tempfile, "colplot" )
r = fwrite( tempfile, "axdelta %axdeltaX" )
r = fwrite( tempfile, "axpos %axposX" )
r = fwrite( tempfile, "axformat %axformX" )
r = fwrite( tempfile, "axis bp" )
r = fwrite( tempfile, "axtitle %axtitleX" )
r = fwrite( tempfile, "axdelta %axdeltaY" )
r = fwrite( tempfile, "axpos %axposY" )
r = fwrite( tempfile, "axformat %axformY" )
r = fwrite( tempfile, "axis lp" )
if (pn = 5) then
r = fwrite( tempfile, "axtitle %axtitleY" )
cif
r = fwrite( tempfile, "axis t" )
r = fwrite( tempfile, "axis r" )
r = fwrite( tempfile, "mm" )
r = fwrite( tempfile, "move %x %y+%Ysizeplot+%titlegap" )
r = fwrite( tempfile, "charheight %chplottitle" )
r = fwrite( tempfile, "text %plottitle(pn)" )
r = fwrite( tempfile, "charheight %chlabels" )
r = fwrite( tempfile, "move %x+%beamDx*%Xsizeplot %y+%beamDy*%Ysizeplot" )
r = fwrite( tempfile, "world" )
r = fwrite( tempfile, "beam %bmmin ARCSEC %bmmaj ARCSEC" )
if (pn = 5) then
r = fwrite( tempfile, "move %raopt %decopt" )
r = fwrite( tempfile, "colour red")
r = fwrite( tempfile, "beam %ellmaj ARCMIN %ellmin ARCMIN %ellpa 0" )
r = fwrite( tempfile, "colour foreground")
cif
plotmade = true
cif
if (pn = 4) then
!--------------------------------------------------
! Does Ascii file with global profile exists?
! first try to open the file including a PGC number.
! If that is not successful, open it without any
! PGC number.
!--------------------------------------------------
setname(pn) = "flux%ugcnr_s60hanp%pgcnr(i)_.dat"
exist(pn) = existf( setname(pn) )
fctor = 1.0
if (not exist(pn)) then
setname(pn) = "flux%ugcnr_s60han.dat"
exist(pn) = existf( setname(pn) )
fctor = 5.0
if (not exist(pn)) then
write "[%progname] Could not open Ascii file [%setname(pn)]"
cif
cif
if (exist(pn)) then
!--------------------------------------------------
! Keep flexibility by determining x and y ranges
! of the global profile data by COLA itself.
!--------------------------------------------------
datfile = fopen( setname(pn), "r" )
first = true
repeat
r = fread( datfile, line )
morelines = (r >= 0)
if (morelines) then
if (first) then
xdatmax = atof( substr1( line, 1 ) )
xdatmin = xdatmax
ydatmax = atof( substr1( line, 2 ) )
ydatmin = ydatmax
first = false
else
xdat = atof( substr1( line, 1 ) )
ydat = atof( substr1( line, 2 ) )
if (xdat < xdatmin) then
xdatmin = xdat
elsif (xdat > xdatmax) then
xdatmax = xdat
cif
if (ydat < ydatmin) then
ydatmin = ydat
elsif (ydat > ydatmax) then
ydatmax = ydat
cif
cif
cif
until (not morelines)
r = fclose( datfile )
xrange = "%xdatmin %xdatmax"
ydatmin = 1.1*fctor*ydatmin
ydatmax = 1.1*fctor*ydatmax
yrange = "%ydatmin %ydatmax"
write "[%progname] Plot %pn is file [%setname(pn)], title: %plottitle(pn)"
r = fwrite( tempfile, "xcolumn file(%setname(pn),1,1:)" )
r = fwrite( tempfile, "ycolumn file(%setname(pn),2,1:)" )
r = fwrite( tempfile, "eycolumn file(%setname(pn),5,1:)" )
r = fwrite( tempfile, "ycolumn %fctor*ycolumn" )
r = fwrite( tempfile, "eycolumn %fctor*eycolumn" )
r = fwrite( tempfile, "xrange %xrange" )
r = fwrite( tempfile, "yrange %yrange" )
r = fwrite( tempfile, "xsize %Xsizeplot" )
r = fwrite( tempfile, "ysize %Ysizeplot" )
r = fwrite( tempfile, "axis bw" )
r = fwrite( tempfile, "axtitle Radial velocity (km/s)" )
r = fwrite( tempfile, "axis lw" )
r = fwrite( tempfile, "axtitle Flux density (mJy)" )
r = fwrite( tempfile, "axis t" )
r = fwrite( tempfile, "axis r" )
r = fwrite( tempfile, "colour blue" )
r = fwrite( tempfile, "points" )
r = fwrite( tempfile, "connect" )
r = fwrite( tempfile, "colour red" )
r = fwrite( tempfile, "errorbar -y y")
r = fwrite( tempfile, "colour foreground" )
r = fwrite( tempfile, "mm" )
r = fwrite( tempfile, "move %x %y+%Ysizeplot+%titlegap" )
r = fwrite( tempfile, "charheight %chplottitle" )
r = fwrite( tempfile, "text %plottitle(pn)" )
r = fwrite( tempfile, "charheight %chlabels" )
r = fwrite( tempfile, "world" )
r = fwrite( tempfile, "move %vsys %ydatmin" )
r = fwrite( tempfile, "justifi 0.5" )
r = fwrite( tempfile, "charheight %charrow" )
r = fwrite( tempfile, "text \\(2262)" )
r = fwrite( tempfile, "move %vsys+(%xdatmax-%xdatmin)/15.0 %ydatmin+(%ydatmax-%ydatmin)/20.0" )
r = fwrite( tempfile, "charheight %chlabels" )
r = fwrite( tempfile, "text V\\dsys\\u" )
r = fwrite( tempfile, "justifi 0.0" )
plotmade = true
cif
cif
if (pn = 6) then
!--------------------------------------------------
! XV plot
!--------------------------------------------------
setname(pn) = "u%ugcnr_s30xv%pgcnr(i)"
exist(pn) = exists( setname(pn) )
if (not exist(pn)) then
write "[%progname] Could not open set [%setname(pn)]"
cif
write "Plot %pn is set [%setname(pn)], title: %plottitle(pn)"
if (exist(pn)) then
nsubs = getsas( "%setname(pn) dec 0", set, subsets )
exist(pn) = (nsubs > 0)
if (not exist(pn)) then
write "[%progname] Something wrong with set [%setname(pn) dec 0]"
cif
cif
if (exist(pn)) then
r = fwrite( tempfile, "inset %setname(pn) dec 0" )
fromheader = (existt("%setname(pn)", subsets(1), "MBOX") > 0)
boxstr = ""
if (fromheader) then
r = rcolumn( "%setname(pn)", subsets(1), "MBOX", "MINBOX", boxstr, 1, 1 )
fromheader = (len(boxstr) > 0)
if (fromheader) then
read "%Plotid(pn): Retrieve BOX from header? [Y]/N" fromheader
cancel fromheader
cif
cif
if (not fromheader) then
"VIEW INSET=%setname(pn) dec 0 CLIP="
"MINBOX INSET=%setname(pn) dec 0 OUTSET= SQUARE=yes"
r = rcolumn( "%setname(pn)", subsets(1), "MBOX", "MINBOX", boxstr, 1, 1 )
cif
boxstr = substr2(boxstr, 1, len(boxstr) )
r = fwrite( tempfile, "box %boxstr" )
r = fwrite( tempfile, "xsize %Xsizeplot" )
r = fwrite( tempfile, "ysize %Ysizeplot" )
!--------------------------------------------------
! At this point, the user wants to see whether he
! needs to adjust the layout of the axis labels
! by adjusting values of the axis label start
! position and delta.
!--------------------------------------------------
fromheader = (existt("%setname(pn)", subsets(1), "AXES") > 0)
offunits = "arcmin"
velunits = "km/s"
unitsdelta = ""
unitspos = ""
axdeltaX = "1"
axdeltaY = ""
! axposX = ""
axposY = ""
axformX = ""
axformY = ""
axtitleX = "Offset from center (%offunits)"
axtitleY = "Radial velocity (km/s)"
if (fromheader) then
ok = true
read "%Plotid(pn): Retrieve AXIS properties from header?: [Y]/N" ok
cancel ok
if (ok) then
r = rcolumn( setname(pn), subsets(1), "AXES", "OFFUNIT", offunits, 1, 1 )
r = rcolumn( setname(pn), subsets(1), "AXES", "VELUNIT", velunits, 1, 1 )
r = rcolumn( setname(pn), subsets(1), "AXES", "AXDELTX", axdeltaX, 1, 1 )
r = rcolumn( setname(pn), subsets(1), "AXES", "AXDELTY", axdeltaY, 1, 1 )
! r = rcolumn( setname(pn), subsets(1), "AXES", "AXPOSX", axposX, 1, 1 )
r = rcolumn( setname(pn), subsets(1), "AXES", "AXPOSY", axposY, 1, 1 )
r = rcolumn( setname(pn), subsets(1), "AXES", "AXFORMX", axformX, 1, 1 )
r = rcolumn( setname(pn), subsets(1), "AXES", "AXFORMY", axformY, 1, 1 )
r = rcolumn( setname(pn), subsets(1), "AXES", "AXTITX", axtitleX, 1, 1 )
r = rcolumn( setname(pn), subsets(1), "AXES", "AXTITY", axtitleY, 1, 1 )
r = len( axdeltaY )
unitsdelta = ""
if (r > 0) then
unitsdelta = "%velunits"
cif
unitspos = ""
if (r > 0) then
unitspos = "%velunits"
cif
cif
cif
repeat
r = list( 0 )
"GPLOT COMMAND=device x11;inset %setname(pn) dec 0;box %boxstr;autoscale;
axdelta %axdeltaX %offunits; axform %axformX;axis bp;axtitle %axtitleX;
axdelta %axdeltaY %unitsdelta; axpos %axposY %unitspos;axform %axformY;axis lp;axtitle %axtitleY;axis t;axis r;quit"
r = list( 1 )
ok = true
read "Axes ok? [Y]/N" ok
cancel ok
if (not ok) then
defstr = substr2(offunits, 1, len(offunits) )
read "Enter units of offset axis: [%defstr]" offunits
cancel offunits
defstr = substr2(axdeltaX, 1, len(axdeltaX) )
offunits = substr2(offunits, 1, len(offunits) )
read "Give delta (%offunits) for X-axis labeling: [%defstr]" axdeltaX
cancel axdeltaX
! defstr = substr2(axposX, 1, len(axposX) )
! read "Give position of 1th X-axis label: [%defstr]" axposX
! cancel axposX
defstr = substr2(axformX, 1, len(axformX) )
read "Give format of X-axis labels: [%defstr]" axformX
cancel axformX
defstr = substr2(axtitleX, 1, len(axtitleX) )
read "Give title along X-axis: [%defstr]" axtitleX
cancel axtitleX
defstr = substr2(axdeltaY, 1, len(axdeltaY) )
read "Give delta (%velunits) for Y-axis labeling: [%defstr]" axdeltaY
cancel axdeltaY
unitsdelta = ""
r = len( axdeltaY )
if (r > 0) then
unitsdelta = "%velunits"
cif
defstr = substr2(axposY, 1, len(axposY) )
read "Give position (%velunits_) of 1th Y-axis label: [%defstr]" axposY
cancel axposY
r = len( axposY )
unitspos = ""
if (r > 0) then
unitspos = "%velunits"
cif
defstr = substr2(axformY, 1, len(axformY) )
read "Give format of Y-axis labels: [%defstr]" axformY
cancel axformY
defstr = substr2(axtitleY, 1, len(axtitleY) )
read "Give title along Y-axis: [%defstr]" axtitleY
cancel axtitleY
cif
until (ok)
!--------------------------------------------------
! Store these values in the header
!--------------------------------------------------
if (fromheader) then
r = deletet( setname(pn), subsets(1), "AXES" )
cif
r = len( axdeltaX )
r = createc( setname(pn), subsets(1), "AXES", "AXDELTX", "CHAR%r",
"Space between major tickmarks in X", " " )
r = wcolumn( setname(pn), subsets(1), "AXES", "AXDELTX", axdeltaX, 1, 1 )
r = len( offunits )
r = createc( setname(pn), subsets(1), "AXES", "OFFUNIT", "CHAR%r",
"Units of offset axis", " " )
r = wcolumn( setname(pn), subsets(1), "AXES", "OFFUNIT", offunits, 1, 1 )
r = len( velunits )
r = createc( setname(pn), subsets(1), "AXES", "VELUNIT", "CHAR%r",
"Units of velocity axis", " " )
r = wcolumn( setname(pn), subsets(1), "AXES", "VELUNIT", velunits, 1, 1 )
! r = len( axposX )
! r = createc( setname(pn), subsets(1), "AXES", "AXPOSX", "CHAR%r",
! "Position of first (complete) label in X", " " )
! r = wcolumn( setname(pn), subsets(1), "AXES", "AXPOSX", axposX, 1, 1 )
r = len( axformX )
r = createc( setname(pn), subsets(1), "AXES", "AXFORMX", "CHAR%r",
"Format of the labels in X", " " )
r = wcolumn( setname(pn), subsets(1), "AXES", "AXFORMX", axformX, 1, 1 )
r = len( axtitleX )
r = createc( setname(pn), subsets(1), "AXES", "AXTITX", "CHAR%r",
"Title along X-axis", " " )
r = wcolumn( setname(pn), subsets(1), "AXES", "AXTITX", axtitleX, 1, 1 )
r = len( axdeltaY )
r = createc( setname(pn), subsets(1), "AXES", "AXDELTY", "CHAR%r",
"Space between major tickmarks in Y", " " )
r = wcolumn( setname(pn), subsets(1), "AXES", "AXDELTY", axdeltaY, 1, 1 )
r = len( axposY )
r = createc( setname(pn), subsets(1), "AXES", "AXPOSY", "CHAR%r",
"Position of first (complete) label in Y", " " )
r = wcolumn( setname(pn), subsets(1), "AXES", "AXPOSY", axposY, 1, 1 )
r = len( axformY )
r = createc( setname(pn), subsets(1), "AXES", "AXFORMY", "CHAR%r",
"Format of the labels in Y", " " )
r = wcolumn( setname(pn), subsets(1), "AXES", "AXFORMY", axformY, 1, 1 )
r = len( axtitleY )
r = createc( setname(pn), subsets(1), "AXES", "AXTITY", "CHAR%r",
"Title along Y-axis", " " )
r = wcolumn( setname(pn), subsets(1), "AXES", "AXTITY", axtitleY, 1, 1 )
axdeltaX = substr2(axdeltaX, 1, len(axdeltaX) )
axdeltaY = substr2(axdeltaY, 1, len(axdeltaY) )
! axposX = substr2(axposX, 1, len(axposX) )
axposY = substr2(axposY, 1, len(axposY) )
axformX = substr2(axformX, 1, len(axformX) )
axformY = substr2(axformY, 1, len(axformY) )
axtitleX = substr2(axtitleX, 1, len(axtitleX) )
axtitleY = substr2(axtitleY, 1, len(axtitleY) )
offunits = substr2(offunits, 1, len(offunits) )
velunits = substr2(velunits, 1, len(velunits) )
r = rcolumn( setname(pn), subsets(1), "LEVTAB", "LEVELS", levels, 1, 1 )
fromheader = (len(levels) > 0)
if (not fromheader) then
r = list( 1 )
"MNMX INSET=%setname(pn) dec 0; BOX=%boxstr CHANGE=YES"
r = list( 0 )
r = rhead( set, subsets(1), "DATAMAX", datamax )
if (r < 0) then
datamax = 0.0
cif
r = rhead( set, subsets(1), "DATAMIN", datamin )
if (r < 0) then
datamin = 0.0
cif
levstep = (datamax - datamin) / numlevels
levels = "%.2datamin_:%.2datamax_:%.2levstep"
cif
levels = substr2( levels, 1, len(levels) )
read "%Plotid(pn): Enter levels for this plot [%levels]" levels
cancel levels
if (fromheader) then
r = deletet( setname(pn), subsets(1), "LEVTAB" )
cif
r = len( levels )
r = createc( setname(pn), subsets(1), "LEVTAB", "LEVELS", "CHAR%r",
"Levels in mapunits for contours", " " )
r = wcolumn( setname(pn), subsets(1), "LEVTAB", "LEVELS", levels, 1, 1 )
r = fwrite( tempfile, "levels %levels" )
r = fwrite( tempfile, "contours" )
r = fwrite( tempfile, "axdelta %axdeltaX %offunits" )
r = fwrite( tempfile, "axis bp" )
!-------------------------------------------------
! The titles along the axes are not the defaults!
!-------------------------------------------------
r = fwrite( tempfile, "axtitle %axtitleX" )
r = fwrite( tempfile, "axdelta %axdeltaY %unitsdelta" )
r = fwrite( tempfile, "axpos %axposY %unitspos" )
r = fwrite( tempfile, "axis lp" )
r = fwrite( tempfile, "axtitle %axtitleY" )
r = fwrite( tempfile, "axis t" )
r = fwrite( tempfile, "axis r" )
r = fwrite( tempfile, "mm" )
r = fwrite( tempfile, "move %x %y+%Ysizeplot+%titlegap" )
r = fwrite( tempfile, "charheight %chplottitle" )
r = fwrite( tempfile, "text %plottitle(pn)" )
r = fwrite( tempfile, "charheight %chlabels" )
r = fwrite( tempfile, "colour background")
r = fwrite( tempfile, "fastyle 1" ) ! solid
r = fwrite( tempfile, "move %x+2.0 %y+%Ysizeplot+%titlegap-5.0" )
r = fwrite( tempfile, "rectangle 15 3.5" )
r = fwrite( tempfile, "move %x+3.0 %y+%Ysizeplot+%titlegap-4.0" )
r = fwrite( tempfile, "colour foreground")
r = fwrite( tempfile, "text P.A.: %ellpa (deg)" )
r = fwrite( tempfile, "move %x+0.8*%Xsizeplot %y+%beamDy*%Ysizeplot" )
r = fwrite( tempfile, "world" )
! Put cross beam
r = fwrite( tempfile, "lwidth 4" )
r = fwrite( tempfile, "beam %bmmaj arcsec 1.2*%chsep km/s 0 0 0 3" )
r = fwrite( tempfile, "lwidth 1" )
r = rhead( set, subsets(1), "CDELT1", cdelt1 )
if (r < 0) then
cdelt1 = 1.0
cif
yp = substr1(boxstr,3)
! FWHM divided by 2
r = fwrite( tempfile, "move %ellmaj/120.0/%cdelt1 %yp km/s")
r = fwrite( tempfile, "justifi 0.5" )
r = fwrite( tempfile, "charheight %charrow" )
r = fwrite( tempfile, "text \\(2262)" )
r = fwrite( tempfile, "move -1.0*%ellmaj/120.0/%cdelt1 %yp km/s")
r = fwrite( tempfile, "charheight %charrow" )
r = fwrite( tempfile, "text \\(2262)" )
r = fwrite( tempfile, "justifi 0.0" )
plotmade = true
cif
cif
!-------------------------------------------------
! Finally, if plot could not be made,
! draw a dummy frame
!-------------------------------------------------
if (not plotmade) then
r = fwrite( tempfile, "box -10 -10 10 10" )
r = fwrite( tempfile, "xsize %Xsizeplot" )
r = fwrite( tempfile, "ysize %Ysizeplot" )
r = fwrite( tempfile, "axis bw" )
r = fwrite( tempfile, "axtitle X" )
r = fwrite( tempfile, "axis lw" )
r = fwrite( tempfile, "axtitle Y" )
r = fwrite( tempfile, "axis t" )
r = fwrite( tempfile, "axis r" )
cif
cfor
cif
!--------------------------------------------------
! Decide whether to do the actual plotting or
! continue with the second page half.
!--------------------------------------------------
if (pagehalf = 1) then
gplotfile = "GPu%ugcnr_p%pgcnr(i)"
psfile = "PSu%ugcnr_p%pgcnr(i)"
newpage = false
read "New page after first page half? Y/[N]" newpage
cancel newpage
if (newpage) then
gplotfile = "%gplotfile.txt"
psfile = "%psfile_.ps"
cif
cif
if (pagehalf = 2) then
newpage = true
gplotfile = "%gplotfile_u%ugcnr_p%pgcnr(i)_.txt"
psfile = "%psfile_u%ugcnr_p%pgcnr(i)_.ps"
cif
if (newpage) then
!--------------------------------------------------
! There are two page halfs or the user wants one
! page half on the current page.
! For this page, create a GPLOT input file(name)
! execute the commands and write the results to
! PostScript file.
!--------------------------------------------------
read "GPLOT input file: [%gplotfile]" gplotfile
cancel gplotfile
gplotfileexist = existf( gplotfile )
if (gplotfileexist) then
lg = deletef( gplotfile )
cif
!--------------------------------------------------
! Ask the name of the PostScript file that will be
! created. The default contains the UGC number and
! PGC number(s). There is no check on the existence
! of the file.
!--------------------------------------------------
read "PostScript output file: [%psfile]" psfile
cancel psfile
if ( existf(psfile) ) then
lg = deletef( psfile )
cif
grdevice = "PCPSFILE"
read "Enter GPLOT device: [%grdevice]: " grdevice
cancel grdevice
grdevice = "%grdevice_/%psfile"
write "[%progname] ============= OUTPUT ACTIONS ================= "
write "[%progname] The GPLOT input file is called [%gplotfile]"
write "[%progname] The PostScript output file is called [%psfile]"
write "[%progname] The page half is %pagehalf"
!--------------------------------------------------
! Copy the contents of the temp file to the GPLOT
! input file.
!--------------------------------------------------
r = fclose( tempfile )
tempfile = fopen( tempname, "r" )
gpfp = fopen( gplotfile, "w" )
r = fwrite( gpfp, "device %grdevice" ) ! Write device first
r = fwrite( gpfp, "xmargin %xmargin" ) ! Set some globals
r = fwrite( gpfp, "ymargin %ymargin" )
r = fwrite( gpfp, "font roman" )
repeat
r = fread( tempfile, copyline )
morelines = (r >= 0)
if (morelines) then
r = fwrite( gpfp, copyline )
cif
until (not morelines)
r = fclose( tempfile )
ok = deletef( tempname )
r = fclose( gpfp )
!--------------------------------------------------
! Execute the commands in the GPLOT Ascii file
!--------------------------------------------------
r = list( gplotoutput ) ! set output level
"GPLOT COMMAND=input %gplotfile;quit" ! start GPLOT
!--------------------------------------------------
! Delete the temp. sets
!--------------------------------------------------
r = list( 1 )
if (pagehalf = 1) then
"DELETE INSET=TMP*; OK=Y;Y;Y;"
else
"DELETE INSET=TMP*; OK=Y;Y;Y;Y;Y;Y;"
cif
viewps = true
read "View the result PostScript file? [Y]/N" viewps
cancel viewps
if (viewps) then
if ( existf(psfile) ) then
r = system( "%psviewer %psfile \&" )
else
write "No PostScript file available!"
cif
cif
cif ! End of 'if newpage'
pagehalf = pagehalf + 1
if (pagehalf = 3) then ! Save method
pagehalf = 1
cif
if (newpage) then
pagehalf = 1
cif
cfor
until (false)
