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)