SUBROUTINE SPLOT( STATUS )
*+
* Name:
* SPLOT
* Purpose:
* Generates plot of spectra with overlaid fits
* Language:
* Starlink Fortran
* Type of Module:
* ASTERIX task
* Invocation:
* CALL SPLOT( STATUS )
* Arguments:
* STATUS = INTEGER (Given and Returned)
* The global status.
* Description:
* Obtains observed data (with instrument response if present) and fit
* model from the environment, and produces a GRAFIX dataset of data with
* model overlaid. This can be done rigorously in the data space, or (in
* the case where model and data spaces are distict; i.e. instrument folding
* has taken place) with some approximation in the model space, or both,
* according to the plot mode selected:
* Mode D/C........................... data/channel space plot
* Mode P ............................ photon space plot
* Mode DP ........................... both data and model plots
* Mode DPR......... ................. data and model with residuals
* Mode R............................. data space residuals (only)
* If multiple datasets are entered then all are displayed.
* For the spectral fitting case, a redshift may be included. The model is
* then interpreted as applying in the source frame, but the model space
* spectrum displayed is that observed in the earth frame (i.e. flux from
* blue shifted source energy bins is mapped onto the corresponding
* unshifted rest frame bins.
* Note - no special treatment is required for data which has been
* likelihood fitted - it is assumed that the dataset given is
* b/g-subtracted, so plotting is unaffected.
* Usage:
* splot {parameter_usage}
* Environment Parameters:
* INP = CHAR (read)
* Input data (either a single dataset or a file of references).
* Z = REAL (read)
* Redshift of source spectrum
* MODEL = CHAR (read)
* Data object containing model specification
* MODE = LITERAL(R)
* Plot mode
* OUT = CHAR (read)
* Output plot file
* Examples:
* {routine_example_text}
* {routine_example_description}
* Pitfalls:
* {pitfall_description}...
* Notes:
* {routine_notes}...
* Prior Requirements:
* {routine_prior_requirements}...
* Side Effects:
* {routine_side_effects}...
* Algorithm:
* Model space dataset values are generated by an approximate method
* (unavoidably), explained in SPLOT_TRANSFM. Note that the
* transformation is only valid if the model is a good fit.
* Formatting of plotted output is expressed in terms of 'plot zones'
* 'plots' and 'graphs'. A 'plot' is a set of axes with one or more
* 'graphs' plotted on it (i.e. overlaying allows several graphs per plot).
* A 'plot zone' is just the collection of plots for a single input
* dataset. The total graphical output from the program consists of one
* p.z. for each dataset.
* Note that the base axes in each plot are set by the observed data. This
* may mean that some of an overlaid fit will disappear off the graph.
* Graph legends and axis labels etc are propagated from the input datasets
* into the plots produced.
* Accuracy:
* {routine_accuracy}
* Timing:
* {routine_timing}
* Implementation Status:
* {routine_implementation_status}
* External Routines Used:
* {name_of_facility_or_package}:
* {routine_used}...
* Implementation Deficiencies:
* Assumes that either ALL or NO datasets are folded.
* Residual plotting not yet implemented.
* References:
* {task_references}...
* Keywords:
* splot, usage:public
* Copyright:
* Copyright (C) University of Birmingham, 1995
* Authors:
* TJP: Trevor Ponman (University of Birmingham)
* RJV: Bob Vallance (ROSAT, University of Birmingham)
* {enter_new_authors_here}
* History:
* 26 Aug 88: Converted to FIT_MODEL, model space plotting working (TJP)
* 1 Dec 88: Residual plotting implemented + various minor bug fixes (TJP)
* 24 Jan 89: QUALITY copied into residual plot (TJP)
* 17 May 89: Safeguarded against undefined returns from MATH_LINTERP (TJP)
* 13 Jun 89: Converted to ASTERIX88 (RJV)
* 21 Jun 89: Final mods and redshifting (TJP)
* 12 Jul 89: Final final mods to deal with spectral sets (RJV/TJP)
* 17 Nov 89: Warning if no response found in spectral case (TJP)
* 5 Dec 89: Photon space residuals expressed as significances (TJP)
* 18 May 90: Trap missing quality in model space resid case (TJP)
* 6 Dec 90: Doesn't copy input quality if not present (DJA)
* 5 Aug 90: QUALITY copied to residuals plot! (TJP)
* 22 Jun 92: Handle case when photon space model (& its variance) = 0 (TJP)
* 3 Jul 92: Change to FIT_DATINGET & _PREDDAT interfaces (l-fitting) (TJP)
* 8 Jul 92: Legend change in plots (TJP)
* 9 Aug 92: Bug in handling Quality for spectral sets fixed (TJP)
* 25 Jun 87: V0.6-1 Original (BHVAD::TJP)
* 23 Aug 88: V0.6-2 Major reorganisation, FIT_PLOT introduced (TJP)
* 28 Nov 88: V0.6-3 Residual plotting option implemented (TJP)
* 21 Jun 89: V1.0-1 ASTERIX88 release, incorporates redshifting (TJP/RJV)
* 7 Dec 90: V1.4-1 Bug with multiple spectra and model compts fixed (TJP)
* 5 Aug 91: V1.5-1 Quality arrays copied to resids plots, INIT added (TJP)
* 24 Jun 92: V1.6-1 FIT_PLOT adapted to call l-fit routines + bug fix (TJP)
* Feb 93: V1.7-0 FIT_PLOT got rid of - interactive plotting (RJV)
* 4 Mar 93: V1.7-1 Use new SFIT_ subroutines (DJA)
* 6 May 93: V1.7-2 Divide by bin width (RJV)
* 8 Sep 93: V1.7-3 Added SPEC_INIT call and removed reference to
* SPEC_CMN_RZ (DJA)
* Nov 93: V1.7-4 Complete rebuild and new features (RJV)
* 19 Mar 94: V1.7-5 Model component plot (RJV)
* 24 Nov 94 : V1.8-0 Now use USI for user interface (DJA)
* 24 Apr 95 : V1.8-1 Updated input and model data interfaces. Use
* ADI style response matrices (DJA)
* 30 Nov 1995 V2.0-0 (DJA):
* Full ADI port
* {enter_changes_here}
* Bugs:
* {note_any_bugs_here}
*-
* Type Definitions:
IMPLICIT NONE ! No implicit typing
* Global Constants:
INCLUDE 'SAE_PAR' ! Standard SAE constants
INCLUDE 'FIT_PAR'
INCLUDE 'ADI_PAR'
* Structure Definitions:
INCLUDE 'FIT_STRUC'
* Global Variables:
INCLUDE 'SPLOT_CMN'
* {global_name}[dimensions] = {data_type} ({global_access_mode})
* [global_variable_purpose]
* Status:
INTEGER STATUS ! Global status
* Local Constants:
CHARACTER*30 VERSION
PARAMETER ( VERSION = 'SPLOT Version V2.0-0' )
* Local Variables:
RECORD /DATASET/ OBDAT(NDSMAX) ! Observed datasets
RECORD /INSTR_RESP/ INSTR(NDSMAX) ! Instrument responses
RECORD /PREDICTION/ PREDDAT(NDSMAX) ! Data predicted by model
RECORD /MODEL_SPEC/ MODEL ! Model specification
RECORD /MODEL_SPEC/ TMODEL ! Term model spec
CHARACTER*80 MODSPEC ! Model spec (for title)
CHARACTER*80 TITLE(2) ! Text specifying dataset and model
CHARACTER*20 DEV ! Graphics device
CHARACTER*20 COMPSPEC(MAXCOMP) ! spec for individual components
REAL PARAM(NPAMAX) ! Model parameters
REAL LB(NPAMAX) ! Parameter lower bounds
REAL UB(NPAMAX) ! Parameter upper bounds
REAL LE(NPAMAX) ! Lower parameter errors
REAL UE(NPAMAX) ! Upper parameter errors
REAL XMIN,XMAX ! Min/max x for current p.z. (NDC)
REAL YMIN,YMAX ! Min/max y for current p.z. (NDC)
REAL ELOW,EUPP ! energy bounds for model plot
REAL EWID ! channel width for model plot
REAL Z ! Redshift [ Eobs/Esource=1/(1+z) ]
INTEGER IERR ! First VEC_ error
INTEGER IFID ! Input dataset id
INTEGER IG ! Current graph number
INTEGER IP ! Current plot number
INTEGER ITERM ! Loop over terms
INTEGER IUX ! Current x posn (in plot zone)
INTEGER IUY ! Current y posn (in plot zone)
INTEGER LDIM(2),UDIM(2) ! Slice specification
INTEGER MFID ! Model spec dataset id
INTEGER MOBJ ! O/p multi-graph
INTEGER NCHAN ! # energy channels for model plot
INTEGER NERR ! # VEC_ errors
INTEGER N ! Dataset index
INTEGER NDMAX,NFMAX
INTEGER NMAX,NEL
INTEGER NGOOD ! # good data points
INTEGER NNDF ! Total # NDFs (curves)
INTEGER NPAR ! No of parameters
INTEGER NPZ ! # plots per plot zone
INTEGER NGZ ! # graphs per plot zone
INTEGER NSET ! # spectra in set
INTEGER NTERM ! # additive terms
INTEGER NZX,NZY ! x,y layout of plot zones
INTEGER OFID ! Output dataset
INTEGER QPTR ! I/p quality pointer
INTEGER SIGNS(MAXCOMP) ! Sign of each term
INTEGER SSCALE ! Statistic scaling
INTEGER TERMS(MAXCOMP,MAXCOMP) ! Individual terms
LOGICAL ACTIVE ! Graphics device active
LOGICAL CHANGE ! Plot mode has changed
LOGICAL COMPLEX ! Complex plot (>1 per page)
LOGICAL FROZEN(NPAMAX) ! Parameter frozen?
LOGICAL SAVE ! Save data to file?
*.
* Check inherited global status.
IF ( STATUS .NE. SAI__OK ) RETURN
* Version id
CALL MSG_PRNT( VERSION )
* Initialise ASTERIX
CALL AST_INIT()
CALL SPEC_INIT( STATUS )
* Set up genus in MODEL structure
MODEL.GENUS = 'SPEC'
* Set plot mode flags in common
CALL SPLOT_PLOTMODE( 'MODE', CHANGE, COMPLEX, STATUS )
* Model only plot?
IF ( MPLOT ) THEN
CALL USI_GET0R( 'ELOW', ELOW, STATUS )
CALL USI_GET0R( 'EUPP', EUPP, STATUS )
CALL USI_GET0I( 'NCHAN', NCHAN, STATUS )
ELSE
* Get observed data and response
CALL USI_ASSOC( 'INP', 'FileSet|BinDS', 'READ', IFID, STATUS )
CALL FIT_GETDAT( ADI__NULLID, IFID, 'SPEC', 1, .FALSE.,
: .FALSE., NDS, OBDAT,
: NGOOD, SSCALE, PREDDAT, INSTR, STATUS )
NDMAX = OBDAT(1).NDAT
NFMAX = PREDDAT(1).NMDAT
END IF
* Get model specification
CALL USI_ASSOC( 'MODEL', '*', 'READ', MFID, STATUS )
CALL FIT_MODGET( MFID, MODEL, NPAR, PARAM, LB, UB, LE, UE, FROZEN,
: STATUS )
MODSPEC = MODEL.SPEC
* Look for redshift
CALL SFIT_GETZ( Z, STATUS )
IF ( MPLOT ) THEN
* Apply redshift
CALL SFIT_APPRED(Z,1,ELOW,EUPP,STATUS)
EWID = (EUPP-ELOW)/REAL(NCHAN)
ELSE
* Apply red-shift and check data structures
CALL SFIT_PRECHK( NDS, Z, PREDDAT, STATUS )
* Loop through input datasets to find maximum data and predicted data sizes
DO N = 1, NDS
NDMAX = MAX(NDMAX,OBDAT(N).NDAT)
NFMAX = MAX(NFMAX,PREDDAT(N).NMDAT)
END DO
* Get dynamic storage
CALL SPLOT_GET_DYN( NDMAX, NFMAX, STATUS )
END IF
* Set up workspace for model stack
IF ( MPLOT ) THEN
CALL DYN_MAPR(1,NCHAN*MAXSTACK,MODEL.STACKPTR,STATUS)
ELSE
NMAX = MAX(NDMAX,NFMAX)
CALL DYN_MAPR(1,NMAX*MAXSTACK,MODEL.STACKPTR,STATUS)
END IF
IF ( STATUS .NE. SAI__OK ) GOTO 9000
IF ( MPLOT ) THEN
* Break up model into additive components
CALL FIT_MSPECFLAT( MODEL, MAXCOMP, NTERM, TERMS, SIGNS,
: STATUS )
* Dynamic storage
CALL SPLOT_MPLOT_GET_DYN( NCHAN, NTERM, STATUS )
CALL ARR_INIT1R( 0.0, NCHAN, %val(MTDPTR), STATUS )
* Get channel boundaries and centres
CALL SPLOT_MPLOT_EAXIS( NCHAN, ELOW, EUPP, %VAL(ELPTR),
: %VAL(EUPTR), %VAL(MCXPTR), STATUS )
* Get data values for each term
DO ITERM = 1, NTERM
* Calculate values in a term
CALL FIT_MCALCTERM( MODEL, NTERM, TERMS, SIGNS, ITERM,
: PARAM, 1, 1, NCHAN, NCHAN, NCHAN+1,
: %VAL(ELPTR), %VAL(EUPTR),
: %VAL(MODEL.STACKPTR), TMODEL,
: %VAL(MCDPTR(ITERM)), STATUS )
* Store term model spec for labelling
COMPSPEC(ITERM) = TMODEL.SPEC
* Normalise
CALL SPLOT_DIVWID2( NCHAN, EWID, %VAL(MCDPTR(ITERM)), STATUS )
* Add to total model (what about -ve)
CALL SPLOT_MPLOT_ADDCOMP( NCHAN, %VAL(MCDPTR(ITERM)),
: %VAL(MTDPTR), STATUS )
* Next term
END DO
ELSE
* Loop through input datasets getting predicted data values
DO N = 1, NDS
CALL FIT_PREDDAT( 1, NDS, OBDAT, INSTR, PREDDAT, MODEL,
: PARAM, N, %VAL(PREDDAT(N).DPTR), STATUS )
END DO
* If no convolution is involved, plot mode 'P' is unavailable
IF ( PPLOT .AND. .NOT. (PREDDAT(1).CONVOLVE) ) THEN
CALL MSG_PRNT('No convolution available - '//
: 'cannot plot in photon space')
PPLOT = .FALSE.
PRPLOT = .FALSE.
END IF
END IF
IF ( STATUS .NE. SAI__OK ) GOTO 9000
* Set up layout parameters
CALL SPLOT_LAYOUT( NPZ, NGZ, NZX, NZY, STATUS )
* Get number of NDFs to store data
IF ( MPLOT ) THEN
NNDF = NTERM + 1
ELSE
NNDF = NGZ * NDS
END IF
* Interactive or save mode
CALL USI_GET0L( 'SAVE', SAVE, STATUS )
IF ( STATUS .NE. SAI__OK ) GOTO 9000
* Output is multi-graph
CALL ADI_NEW0( 'MultiGraph', MOBJ, STATUS )
* Create file in save mode
IF ( SAVE ) THEN
CALL USI_CREAT( 'OUT', MOBJ, OFID, STATUS )
CALL MSG_PRNT('Creating multiple dataset...')
* Otherwise write to graphics device
ELSE
* Temporary internal file
CALL ADI_TEMP( OFID, STATUS )
CALL ADI_SETLNK( MOBJ, OFID, STATUS )
OFID = MOBJ
* Open device if not already open
CALL GDV_STATUS( ACTIVE, STATUS )
IF ( .NOT. ACTIVE ) THEN
CALL USI_GET0C( 'DEV', DEV, STATUS )
CALL GDV_OPEN( DEV, 1, 1, STATUS )
END IF
END IF
* Create structure for holding NDFs
CALL GMI_CREMULT( OFID, NNDF, STATUS )
* Connect to local graphics control
CALL GCB_LCONNECT( STATUS )
IF ( STATUS .NE. SAI__OK ) GOTO 9000
* Finally do plotting
* model only plot
IF ( MPLOT ) THEN
CALL SPLOT_MPLOT_SETUP( NCHAN, NTERM, MODSPEC, COMPSPEC,
: STATUS )
CALL SPLOT_MPLOT_COPY( OFID, NCHAN, NTERM, COMPSPEC, SAVE,
: STATUS )
ELSE
* Loop through datasets
DO N = 1, NDS
* Current graph number
IG = 1 + (N-1)*NGZ
* Current plot number
IF ( OVERLAY ) THEN
IP = 1
ELSE
IP = 1+(N-1)*NPZ
END IF
* Current X and Y position in plot zones
IF ( OVERLAY ) THEN
IUX = 1
IUY = 1
ELSE
IUX = 1 + MOD(N-1,NZX)
IUY = 1 + INT((N-1)/NZX)
END IF
* Window for current plot zone (in Normalised Device Coords)
XMIN = REAL(IUX-1)/NZX
XMAX = REAL(IUX)/NZX
YMIN = REAL(IUY-1)/NZY
YMAX = REAL(IUY)/NZY
* Number of values in model space, and address of model space data
NFIT = PREDDAT(N).NMDAT
IFPTR = PREDDAT(N).DPTR
* Number of values in channel space, and address of channel space data
DSID = OBDAT(N).D_ID
NVAL = OBDAT(N).NDAT
IDPTR = OBDAT(N).DPTR
* Number of spectral sets in data
NSET = OBDAT(N).SETINDEX
* Quality and variance present in fit data?
QOK = ((OBDAT(N).QPTR).NE.0)
VOK = ((OBDAT(N).VPTR).NE.0)
* Spectral set?
IF ( NSET .GT. 0 ) THEN
LDIM(1)=1
LDIM(2)=NSET
UDIM(1)=NVAL
UDIM(2)=NSET
NEL = NSET * NVAL
ELSE ! single dataset
LDIM(1)=1
UDIM(1)=NVAL
NEL = NVAL
END IF
* Make copy of data and model
CALL ARR_COP1R( NVAL, %VAL(IDPTR), %VAL(DDPTR), STATUS )
CALL ARR_COP1R( NVAL, %VAL(IFPTR), %VAL(DFDPTR), STATUS )
* Copy data variance if present
IF ( VOK ) THEN
IVPTR = OBDAT(N).VPTR
CALL ARR_COP1R( NVAL, %VAL(IVPTR), %VAL(DVPTR), STATUS )
* If not present, use input data assuming Poisson statistics
ELSE
CALL DYN_MAPR( 1, NEL, DVPTR, STATUS )
CALL VEC_ABSR( .FALSE., NEL, %VAL(IDPTR), %VAL(DVPTR),
: IERR, NERR, STATUS )
END IF
* Copy data quality if present
IF ( QOK ) THEN
CALL BDI_MAPUB( DSID, 'Quality', 'READ', QPTR, STATUS )
CALL DYN_MAPB( 1, NEL, DQPTR, STATUS )
CALL ARR_SLCOPB( OBDAT(N).NDIM, OBDAT(N).IDIM, %VAL(QPTR),
: LDIM, UDIM, %VAL(DQPTR), STATUS )
CALL BDI_GET0UB( DSID, 'QualityMask', MASK, STATUS )
ELSE
DQPTR = 0
END IF
* Access data channel axis values
CALL BDI_AXMAPR( DSID, 1, 'Data', 'READ', DXPTR, STATUS )
CALL BDI_AXMAPR( DSID, 1, 'Width', 'READ', DWPTR, STATUS )
* Set up data for plots
CALL SPLOT_GEN_SETUP( INSTR(N), PREDDAT(N), Z, STATUS )
* Construct plot title
CALL SPLOT_TITLE( OBDAT(N).DATNAME, MODSPEC, TITLE )
* Progress report...
IF ( SAVE ) THEN
CALL MSG_SETI( 'DS', N )
CALL MSG_PRNT( 'Dataset ^DS' )
END IF
* Data space plot
IF ( DPLOT ) THEN
CALL SPLOT_DPLOT_SETUP( N, NPZ, XMIN, XMAX, YMIN, YMAX,
: TITLE, STATUS )
CALL SPLOT_DPLOT_COPY( OFID, IP, IG, N, SAVE, STATUS )
* Update plot and graph nos.
IP = IP + 1
IG = IG + 2
END IF
* Data space residuals
IF ( DRPLOT ) THEN
CALL SPLOT_DRPLOT_SETUP( N, NPZ, XMIN, XMAX, YMIN, YMAX,
: STATUS )
CALL SPLOT_DRPLOT_COPY( OFID, IP, IG, N, SAVE, STATUS )
* Update plot and graph nos.
IP = IP + 1
IG = IG + 1
END IF
* Photon space plot
IF ( PPLOT ) THEN
CALL SPLOT_PPLOT_SETUP( N, NPZ, XMIN, XMAX, YMIN, YMAX,
: TITLE, STATUS )
CALL SPLOT_PPLOT_COPY( OFID, IP, IG, N, SAVE, STATUS )
* Update plot and graph nos.
IP = IP + 1
IG = IG + 2
END IF
* Photon space residuals
IF ( PRPLOT ) THEN
CALL SPLOT_PRPLOT_SETUP( N, NPZ, XMIN, XMAX, YMIN, YMAX,
: STATUS )
CALL SPLOT_PRPLOT_COPY( OFID, IP, IG, N, SAVE, STATUS )
* Update plot and graph nos.
IP = IP + 1
IG = IG + 1
END IF
* End of dataset loop
END DO
* Free workspace
CALL SPLOT_RELEASE_DYN( STATUS )
END IF
* Release model stack
CALL DYN_UNMAP( MODEL.STACKPTR, STATUS )
* Release output
IF ( SAVE ) THEN
CALL USI_ANNUL( 'OUT', STATUS )
* If in interactive mode - plot
ELSE
CALL GDRAW_SUB( OFID, STATUS )
END IF
* Exit
9000 IF ( .NOT. ACTIVE ) THEN
CALL GDV_CLOSE( STATUS )
END IF
CALL GCB_DETACH( STATUS )
CALL AST_CLOSE()
END
*+
SUBROUTINE SPLOT_GET_DYN(ND,NF,STATUS)
* Description :
* Type definitions :
IMPLICIT NONE
* Import :
INTEGER ND,NF
* Import-Export :
* Export :
* Status :
INTEGER STATUS
* Global constants :
INCLUDE 'SAE_PAR'
INCLUDE 'FIT_PAR'
* Global variables :
INCLUDE 'SPLOT_CMN'
* Local variables :
*-
* Check inherited global status.
IF ( STATUS .NE. SAI__OK ) RETURN
CALL DYN_MAPR(1,ND,DDPTR,STATUS)
CALL DYN_MAPR(1,ND,DVPTR,STATUS)
CALL DYN_MAPR(1,ND,DRDPTR,STATUS)
CALL DYN_MAPR(1,ND,DFDPTR,STATUS)
CALL DYN_MAPR(1,ND,PDPTR,STATUS)
CALL DYN_MAPR(1,ND,PVPTR,STATUS)
CALL DYN_MAPR(1,ND,PXPTR,STATUS)
CALL DYN_MAPR(1,ND,PWPTR,STATUS)
CALL DYN_MAPR(1,NF,PFDPTR,STATUS)
CALL DYN_MAPR(1,NF,PFXPTR,STATUS)
CALL DYN_MAPR(1,ND,PRDPTR,STATUS)
CALL DYN_MAPR(1,ND,PRVPTR,STATUS)
CALL DYN_MAPB(1,ND,PRQPTR,STATUS)
CALL DYN_MAPR(1,ND,PIDPTR,STATUS)
IF (STATUS.NE.SAI__OK) THEN
CALL AST_REXIT('SPLOT_GET_DYN',STATUS)
END IF
END
*+
SUBROUTINE SPLOT_RELEASE_DYN(STATUS)
* Description :
* Type definitions :
IMPLICIT NONE
* Import :
* Import-Export :
* Export :
* Status :
INTEGER STATUS
* Global constants :
INCLUDE 'SAE_PAR'
INCLUDE 'FIT_PAR'
* Global variables :
INCLUDE 'SPLOT_CMN'
* Local variables :
*-
CALL DYN_UNMAP(DDPTR,STATUS)
IF ( .NOT. VOK ) THEN
CALL DYN_UNMAP(DVPTR,STATUS)
END IF
CALL DYN_UNMAP(DRDPTR,STATUS)
CALL DYN_UNMAP(DFDPTR,STATUS)
CALL DYN_UNMAP(PDPTR,STATUS)
CALL DYN_UNMAP(PVPTR,STATUS)
CALL DYN_UNMAP(PXPTR,STATUS)
CALL DYN_UNMAP(PWPTR,STATUS)
CALL DYN_UNMAP(PFDPTR,STATUS)
CALL DYN_UNMAP(PFXPTR,STATUS)
CALL DYN_UNMAP(PRDPTR,STATUS)
CALL DYN_UNMAP(PRVPTR,STATUS)
CALL DYN_UNMAP(PRQPTR,STATUS)
CALL DYN_UNMAP(PIDPTR,STATUS)
END
*+
SUBROUTINE SPLOT_MPLOT_GET_DYN(NCHAN,NTERM,STATUS)
* Description :
* Type definitions :
IMPLICIT NONE
* Import :
INTEGER NCHAN,NTERM
* Import-Export :
* Export :
* Status :
INTEGER STATUS
* Global constants :
INCLUDE 'SAE_PAR'
INCLUDE 'FIT_PAR'
* Global variables :
INCLUDE 'SPLOT_CMN'
* Local variables :
INTEGER I
*-
* Check inherited global status.
IF ( STATUS .NE. SAI__OK ) RETURN
CALL DYN_MAPR(1,NCHAN+1,ELPTR,STATUS)
CALL DYN_MAPR(1,NCHAN+1,EUPTR,STATUS)
CALL DYN_MAPR(1,NCHAN,MCXPTR,STATUS)
DO I = 1, NTERM
CALL DYN_MAPR(1,NCHAN,MCDPTR(I),STATUS)
END DO
CALL DYN_MAPR(1,NCHAN,MTDPTR,STATUS)
IF ( STATUS .NE. SAI__OK ) THEN
CALL AST_REXIT('SPLOT_MPLOT_GET_DYN',STATUS)
END IF
END
*+
SUBROUTINE SPLOT_MPLOT_RELEASE_DYN(NTERM,STATUS)
* Description :
* Type definitions :
IMPLICIT NONE
* Import :
INTEGER NTERM
* Import-Export :
* Export :
* Status :
INTEGER STATUS
* Global constants :
INCLUDE 'SAE_PAR'
INCLUDE 'FIT_PAR'
* Global variables :
INCLUDE 'SPLOT_CMN'
* Local variables :
INTEGER I
*-
CALL DYN_UNMAP(ELPTR,STATUS)
CALL DYN_UNMAP(EUPTR,STATUS)
CALL DYN_UNMAP(MCXPTR,STATUS)
DO I=1,NTERM
CALL DYN_UNMAP(MCDPTR(I),STATUS)
END DO
CALL DYN_UNMAP(MTDPTR,STATUS)
END
*+ SPLOT_MODRES - calc. model space residuals
SUBROUTINE SPLOT_MODRES(N,D,FIT,V,DR,VR,QR,STATUS)
* Description :
* Type definitions :
IMPLICIT NONE
* Import :
INTEGER N ! No of data values
REAL D(*) ! Data array
REAL FIT(*)
REAL V(*)
* Import-Export :
* Export :
REAL DR(*)
REAL VR(*)
BYTE QR(*)
* Status :
INTEGER STATUS
* Functions :
BYTE BIT_ORUB
* Local constants :
INCLUDE 'QUAL_PAR'
INCLUDE 'SAE_PAR'
* Local variables :
INTEGER I
LOGICAL WARN
*-
* Check inherited global status.
IF ( STATUS .NE. SAI__OK ) RETURN
WARN=.FALSE.
DO I=1,N
IF (V(I).GT.0.0) THEN
DR(I)=(D(I)-FIT(I))/SQRT(V(I))
ELSE
DR(I) = 0.0
QR(I) = BIT_ORUB(QR(I),QUAL__ARITH)
WARN=.TRUE.
END IF
VR(I)=1.0
END DO
IF ( WARN ) THEN
CALL MSG_PRNT(' ')
CALL MSG_PRNT('Warning - photon space model is zero for '//
: 'some channel energies')
CALL MSG_PRNT(' ')
END IF
END
*+ SPLOT_RESID - calc. data space residuals
SUBROUTINE SPLOT_RESID(N,D,FIT,DR,STATUS)
* Description :
* Type definitions :
IMPLICIT NONE
* Import :
INTEGER N ! No of data values
REAL D(*) ! Data array
REAL FIT(*)
* Import-Export :
* Export :
REAL DR(*)
* Status :
INTEGER STATUS
* Local constants :
INCLUDE 'SAE_PAR'
* Local variables :
INTEGER I
*-
IF (STATUS.EQ.SAI__OK) THEN
DO I=1,N
DR(I)=D(I) - FIT(I)
END DO
END IF
END
*+ SPLOT_CHEN - Derives channel centres & widths from bounds
SUBROUTINE SPLOT_CHEN(NDAT,CBOUND,CHEN,CHWID,STATUS)
* Description :
* History :
* 24 Aug 88: Original (TJP)
* 21 Jun 89: ASTERIX88; half-widths -> widths (RJV)
* Type definitions :
IMPLICIT NONE
* Import :
INTEGER NDAT ! No of energy space channels
REAL CBOUND(0:NDAT) ! Channel bound energies
* Import-Export :
* Export :
REAL CHEN(NDAT) ! Channel centre energies
REAL CHWID(NDAT) ! Channel widths
* Status :
INTEGER STATUS
* Global constants :
INCLUDE 'SAE_PAR'
* Local variables :
INTEGER I
*-
IF (STATUS.EQ.SAI__OK) THEN
DO I=1,NDAT
CHEN(I)=(CBOUND(I-1)+CBOUND(I))/2
CHWID(I)=CBOUND(I)-CBOUND(I-1)
END DO
END IF
END
*+ SPLOT_NORM - Normalises flux array w.r.t. energy
SUBROUTINE SPLOT_NORM(NDAT,FLUX,LB,UB,Z,SPECFLUX,STATUS)
* Description :
* Normalises flux to spectral flux observed at earth. Since the
* energy bin bounds are those applying at the source the interval
* between them is decreased by a factor 1/(1+z) at earth, and the
* spectral flux (photon/(cm**2*s*keV)) consequently increases by a
* factor 1+z.
* History :
* 18 Aug 88: Original (TJP)
* 21 Jun 89: Redshifting added (TJP)
* Type definitions :
IMPLICIT NONE
* Import :
INTEGER NDAT ! No of energy space channels
REAL FLUX(NDAT) ! Flux in each channel
REAL LB(NDAT) ! Lower channel bounds
REAL UB(NDAT) ! Upper channel bounds
REAL Z ! Redshift [ Eobs/Esource=1/(1+z) ]
* Import-Export :
* Export :
REAL SPECFLUX(NDAT) ! Flux per unit energy in each channel
* Status :
INTEGER STATUS
* Global constants :
INCLUDE 'SAE_PAR'
* Local variables :
INTEGER I
*-
IF (STATUS.EQ.SAI__OK) THEN
DO I=1,NDAT
SPECFLUX(I)=(1+Z)*FLUX(I)/(UB(I)-LB(I))
END DO
END IF
END
*+ SPLOT_TRANSFM - Transforms observed data into approx model space values
SUBROUTINE SPLOT_TRANSFM(NDAT,OBDAT,PREDDAT,ERR,OBVAR,
: PREDMOD,TRANDAT,TRANVAR,STATUS)
* Description :
* The observed data values OBDAT are transformed approximately back into
* the model space (e.g. channel spectrum transformed to photon spectrum)
* using the prescription:
* TRANDAT = OBDAT*PREDMOD/PREDDAT
* i.e. the returned values are the model space values from the FIT model
* (interpolated to give values at the data channel energies) scaled by
* the ratio of the observed to predicted intensities in each data channel.
* The PREDMOD values should already be normalised to unit model bin width.
* This results, in the spectral case, in output values in
* photon/(cm**2*sec*keV), rather than photon/(cm**2*sec).
* Where errors are present they are scaled similarly, so that
* TRANVAR=OBVAR*(TRANDAT/OBDAT)**2
* (Note that errors are represented as variances in ASTERIX88)
* It is assumed that predicted data values will be positive, and to
* avoid floating zero divides, all are set to be >=PREDLIM.
* History :
* 25 Jun 87: Original (TJP)
* 25 Nov 88: Proofed against PREDDAT=0 (TJP)
* 21 Jun 89: ASTERIX88 version, errors->variances
* Type definitions :
IMPLICIT NONE
* Import :
INTEGER NDAT ! No of spectral channels
REAL OBDAT(*) ! Observed spectrum
REAL PREDDAT(*) ! Predicted spectrum
LOGICAL ERR ! Errors available?
REAL OBVAR(*) ! Data variances
REAL PREDMOD(*) ! Model space prediction
* Import-Export :
* Export :
REAL TRANDAT(*) ! Transformed data values (i.e. in model space)
REAL TRANVAR(*) ! Transformed variances
* Status :
INTEGER STATUS
* Global Constants :
INCLUDE 'SAE_PAR'
* Local Constants :
REAL PREDLIM ! Minimum value allowed for predicted data
PARAMETER (PREDLIM=1.0E-10)
* Local variables :
INTEGER I
REAL PDAT ! PREDDAT value limited to be >=PREDLIM
*-
IF (STATUS.EQ.SAI__OK) THEN
DO I=1,NDAT
IF(PREDDAT(I).GE.PREDLIM)THEN
PDAT=PREDDAT(I)
ELSE
PDAT=PREDLIM
END IF
TRANDAT(I)=OBDAT(I)*PREDMOD(I)/PDAT
IF(ERR)THEN
TRANVAR(I)=OBVAR(I)*(PREDMOD(I)/PDAT)**2
END IF
END DO
END IF
END
*+ SPLOT_ZERO - Zeroes any undefined values in model space array
SUBROUTINE SPLOT_ZERO(N,ARRAY,STATUS)
* Description :
* MATH_LINTERP returns values VAL__BADR if asked for a value outside
* its interpolation range. This can occur if some channel energies lie
* outside the energy range in the response matrix (which they shouldn't
* really...). This routine simply zeroes these on the basis that they
* ought not to contribute to the source counts. This will result in
* the corresponding transformed data and errors being zeroed too.
* History :
* 17 May 89: Original (TJP)
* Type definitions :
IMPLICIT NONE
* Import :
INTEGER N ! No of array elements
* Import-Export :
REAL ARRAY(N) ! Array
* Export :
* Status :
INTEGER STATUS
* Global constants :
INCLUDE 'SAE_PAR'
INCLUDE 'PRM_PAR'
* Local variables :
INTEGER I
*-
IF (STATUS.EQ.SAI__OK) THEN
DO I=1,N
IF(ARRAY(I).LE.VAL__BADR)THEN
ARRAY(I)=0.0
END IF
END DO
END IF
END
SUBROUTINE SPLOT_DIVWID(N,W,D,STATUS)
IMPLICIT NONE
INCLUDE 'SAE_PAR'
INTEGER STATUS
INTEGER N
REAL W(*)
REAL D(*)
INTEGER I
IF (STATUS.EQ.SAI__OK) THEN
DO I=1,N
IF (W(I).GT.0.0) THEN
D(I)=D(I)/W(I)
END IF
END DO
END IF
END
SUBROUTINE SPLOT_DIVWID2(N,W,D,STATUS)
IMPLICIT NONE
INCLUDE 'SAE_PAR'
INTEGER STATUS
INTEGER N
REAL W
REAL D(*)
INTEGER I
IF (STATUS.EQ.SAI__OK.AND.W.GT.0.0) THEN
DO I=1,N
D(I)=D(I)/W
END DO
END IF
END
*+
SUBROUTINE SPLOT_PLOTMODE(PAR,CHANGE,COMPLEX,STATUS)
* Description :
* D = plot observed data (ie counts)
* C = same as D
* P = plot photons (ie convolved data - if convolution available)
* R = plot residuals
* E = all x-axes to be energy (ie not channel number)
* L = all axes log scale - by default photon plot axes are log
* M = plot model by itself, showing each component
* O = overlay multiple datasets, instead of separate plots
* Type definitions :
IMPLICIT NONE
* Global constants :
INCLUDE 'SAE_PAR'
INCLUDE 'FIT_PAR'
* Global variables :
INCLUDE 'SPLOT_CMN'
* Structure definitions :
* Import :
CHARACTER*(*) PAR
* Import-Export :
* Export :
LOGICAL CHANGE
LOGICAL COMPLEX
* Status :
INTEGER STATUS
* Function declarations :
* Local constants :
* Local variables :
CHARACTER*10 MODE,LAST
SAVE LAST
INTEGER NP
LOGICAL D,P,R,E,M,L,O
*-
IF (STATUS.EQ.SAI__OK) THEN
* get plotting mode
MODE=' '
DO WHILE (MODE.EQ.' '.AND.STATUS.EQ.SAI__OK)
CALL USI_DEF0C(PAR,LAST,STATUS)
CALL USI_GET0C(PAR,MODE,STATUS)
CALL CHR_UCASE(MODE)
IF (MODE(:1).EQ.'H'.OR.MODE.EQ.' '.AND.
: STATUS.EQ.SAI__OK) THEN
CALL SPLOT_PLOTMODE_HELP()
MODE=' '
CALL USI_CANCL(PAR,STATUS)
END IF
END DO
IF (STATUS.EQ.SAI__OK) THEN
* see if mode has changed
IF (MODE.NE.LAST) THEN
CHANGE=.TRUE.
LAST=MODE
END IF
DPLOT=.FALSE.
PPLOT=.FALSE.
DRPLOT=.FALSE.
PRPLOT=.FALSE.
MPLOT=.FALSE.
EAXIS=.FALSE.
LOGAX=.FALSE.
OVERLAY=.FALSE.
E=(INDEX(MODE,'E').NE.0)
R=(INDEX(MODE,'R').NE.0)
M=(INDEX(MODE,'M').NE.0)
D=(INDEX(MODE,'D').NE.0.OR.INDEX(MODE,'C').NE.0)
P=(INDEX(MODE,'P').NE.0)
O=(INDEX(MODE,'O').NE.0)
L=(INDEX(MODE,'L').NE.0)
EAXIS=E
NP=0
IF (MODE.EQ.' ') THEN
CALL MSG_PRNT('No plot mode specified - '//
: 'defaulting to data plot')
DPLOT=.TRUE.
NP=NP+1
ELSEIF (M) THEN
MPLOT=.TRUE.
IF (D.OR.P.OR.R) THEN
CALL MSG_PRNT('Cannot mix model plot with others')
CALL MSG_PRNT(' - other modes will be ignored')
END IF
NP=NP+1
ELSE
IF (R) THEN
IF (D) THEN
DRPLOT=.TRUE.
NP=NP+1
END IF
IF (P) THEN
PRPLOT=.TRUE.
NP=NP+1
END IF
IF (.NOT.(D.OR.P)) THEN
DRPLOT=.TRUE.
NP=NP+1
END IF
END IF
IF (D) THEN
DPLOT=.TRUE.
NP=NP+1
END IF
IF (P) THEN
PPLOT=.TRUE.
NP=NP+1
END IF
IF (O) THEN
OVERLAY=.TRUE.
EAXIS=.TRUE.
END IF
END IF
IF (L) THEN
LOGAX=.TRUE.
END IF
COMPLEX=(NP.NE.1.AND..NOT.OVERLAY)
END IF
END IF
END
*+
SUBROUTINE SPLOT_PLOTMODE_HELP()
INTEGER NLINE
PARAMETER (NLINE=7)
INTEGER I
CHARACTER*79 TEXT(NLINE)
DATA TEXT/
: ' D (or C) - plot observed data (counts)',
: ' P - plot convolved data (photons)',
: ' R - plot residuals',
: ' E - all x-axes plotted as energy',
: ' L - all axes plotted with log scale',
: ' O - (multiple data overlayed on same axes)-soon',
: ' M - plot model only, showing components'/
*-
CALL MSG_BLNK()
DO I=1,NLINE
CALL MSG_PRNT(TEXT(I))
END DO
CALL MSG_BLNK()
END
*+
SUBROUTINE SPLOT_TITLE(DATNAME,MODSPEC,TITLE)
* Description :
* Type definitions :
IMPLICIT NONE
* Global constants :
* Global variables :
* Structure definitions :
* Import :
CHARACTER*(*) DATNAME,MODSPEC
* Import-Export :
* Export :
CHARACTER*(*) TITLE(2)
* Status :
* Function declarations :
INTEGER CHR_LEN
* Local constants :
* Local variables :
INTEGER I,L
*-
L=CHR_LEN(DATNAME)
I=L
* strip off directory spec.
DO WHILE (.NOT.(DATNAME(I:I).EQ.']'.OR.DATNAME(I:I).EQ.'/')
: .AND.I.GT.1)
I=I-1
END DO
IF (I.GT.1) THEN
I=I+1
END IF
TITLE(1)='Dataset: '//DATNAME(I:)
TITLE(2)='Model: '//MODSPEC
END
*+
SUBROUTINE SPLOT_LAYOUT(NPZ,NGZ,NZX,NZY,STATUS)
* Description :
* Type definitions :
IMPLICIT NONE
* Global constants :
INCLUDE 'SAE_PAR'
INCLUDE 'FIT_PAR'
* Global variables :
INCLUDE 'SPLOT_CMN'
* Structure definitions :
* Import :
* Import-Export :
* Export :
INTEGER NPZ,NGZ,NZX,NZY
* Status :
INTEGER STATUS
* Function declarations :
* Local constants :
* Local variables :
INTEGER NXZ,NYZ
*-
IF (STATUS.EQ.SAI__OK) THEN
* just plotting model
IF (MPLOT) THEN
NPZ=1
NGZ=1
NXZ=1
NYZ=1
NZX=1
NZY=1
ELSE
* Number of plots per plot zone
NPZ=0
NGZ=0
IF (DPLOT) THEN
NPZ=NPZ+1 ! one plot
NGZ=NGZ+2 ! two curves
END IF
IF (PPLOT) THEN
NPZ=NPZ+1 ! same as DPLOT
NGZ=NGZ+2
END IF
IF (DRPLOT) THEN
NPZ=NPZ+1 ! one curve per plot
NGZ=NGZ+1
END IF
IF (PRPLOT) THEN
NPZ=NPZ+1
NGZ=NGZ+1
END IF
* Layout within each plot zone (dependent on plot mode)
IF(NPZ.EQ.1)THEN
NXZ=1
NYZ=1
ELSEIF(NPZ.EQ.2)THEN
NXZ=1
NYZ=2
ELSEIF(NPZ.EQ.4)THEN
NXZ=2
NYZ=2
END IF
* Layout of plot zones - one per dataset unless overlayed
NZX=1
NZY=1
IF (.NOT.OVERLAY) THEN
DO WHILE(NZX*NZY.LT.NDS)
* nx no more than ny+2
IF(NZX.GE.NZY+2)THEN
NZY=NZY+1
ELSE
NZX=NZX+1
END IF
END DO
END IF
END IF
END IF
END
*+
SUBROUTINE SPLOT_GEN_SETUP(INSTR,PREDDAT,Z,STATUS)
* Description :
* Type definitions :
IMPLICIT NONE
* Global constants :
INCLUDE 'SAE_PAR'
INCLUDE 'FIT_PAR'
* Global variables :
INCLUDE 'SPLOT_CMN'
* Structure definitions :
INCLUDE 'FIT_STRUC'
* Import :
RECORD /INSTR_RESP/ INSTR
RECORD /PREDICTION/ PREDDAT
REAL Z
* Import-Export :
* Export :
* Status :
INTEGER STATUS
* Function declarations :
* Local constants :
REAL YDEC ! No of decades allowed in range of
PARAMETER (YDEC=3.1) ! data
* Local variables :
INTEGER CBPTR,CIPTR,EIPTR,RPTR
INTEGER NACT,NRESP
LOGICAL CBOUNDS
*-
* Check inherited global status.
IF ( STATUS .NE. SAI__OK ) RETURN
* Divide data and fit by bin width to normalise
CALL SPLOT_DIVWID(NVAL,%val(DWPTR),%val(DDPTR),STATUS)
CALL SPLOT_DIVWID(NVAL,%val(DWPTR),%val(DVPTR),STATUS)
CALL SPLOT_DIVWID(NVAL,%val(DWPTR),%val(DVPTR),STATUS)
CALL SPLOT_DIVWID(NVAL,%val(DWPTR),%val(DFDPTR),STATUS)
* Get fit data
CALL ADI_CGET1R( INSTR.R_ID, 'Energy', NFIT, %val(PFXPTR),
: NACT, STATUS )
CALL SPLOT_NORM( NFIT, %val(PREDDAT.MPTR),
: %val(PREDDAT.MLBNDPTR),%val(PREDDAT.MUBNDPTR),Z,
: %val(PFDPTR),STATUS)
* determine channel energies & widths
* look for channel bound information
CALL ADI_CSIZE( INSTR.R_ID, 'Channels', NACT, STATUS )
CALL ADI_CMAPR( INSTR.R_ID, 'Channels', 'READ', CBPTR, STATUS )
IF(STATUS.NE.SAI__OK)THEN
CALL ERR_ANNUL(STATUS)
CBOUNDS=.FALSE.
ELSE IF(NACT.NE.NVAL+1)THEN
CALL MSG_PRNT('Channel bounds array is incorrect length'
: //' and will be ignored')
CBOUNDS=.FALSE.
ELSE
CBOUNDS=.TRUE.
* derive centres and widths from bounds
CALL SPLOT_CHEN(NVAL,%val(CBPTR),%val(PXPTR),
: %val(PWPTR),STATUS)
CALL ADI_CUNMAP( INSTR.R_ID, 'Channels', CBPTR, STATUS )
END IF
* No bounds available, so derive them by interpoln from response matrix.
* Can only do this for ASTERIX style responses - not likely to be a
* problem for non-ASTERIX responses anyhow.
IF ( .NOT. CBOUNDS ) THEN
CALL ADI_CSIZE( INSTR.R_ID, 'RMF', NRESP, STATUS )
CALL ADI_CMAPI( INSTR.R_ID, 'ChannelIndices', 'READ', CIPTR,
: STATUS )
CALL ADI_CMAPI( INSTR.R_ID, 'EnergyIndices', 'READ', EIPTR,
: STATUS )
CALL ADI_CMAPR( INSTR.R_ID, 'RMF', 'READ', RPTR, STATUS )
CALL MSG_PRNT('Deriving channel energies & widths '
: //'from response matrix')
CALL SPEC_CHEN( NVAL, NRESP, %VAL(CIPTR), %VAL(EIPTR),
: %VAL(RPTR), %val(PFXPTR), %val(PXPTR),
: %val(PWPTR), STATUS )
CALL ADI_CUNMAP( INSTR.R_ID, 'ChannelIndices', CIPTR, STATUS )
CALL ADI_CUNMAP( INSTR.R_ID, 'EnergyIndices', EIPTR, STATUS )
CALL ADI_CUNMAP( INSTR.R_ID, 'RMF', RPTR, STATUS )
END IF
* Interpolate to give predicted model space fluxes at channel energies
CALL MATH_LINTERP(NFIT,%val(PFXPTR),%val(PFDPTR),
: NVAL,%val(PXPTR),%val(PIDPTR),STATUS)
* Safeguard against bad (undefined) model space values (from MATH_LINTERP)
CALL SPLOT_ZERO(NVAL,%val(PIDPTR),STATUS)
* Model space data (i.e. transformed back using fitted model)
CALL SPLOT_TRANSFM(NVAL,%val(IDPTR),%val(IFPTR),
: .TRUE.,%val(IVPTR),%val(PIDPTR),
: %val(PDPTR),%val(PVPTR),STATUS)
END
*+
SUBROUTINE SPLOT_DPLOT_SETUP(IDS,NPZ,XMIN,XMAX,YMIN,YMAX,TITLE,
: STATUS)
* Description :
* Type definitions :
IMPLICIT NONE
* Global constants :
INCLUDE 'SAE_PAR'
INCLUDE 'FIT_PAR'
* Global variables :
* Structure definitions :
* Import :
INTEGER IDS,NPZ
REAL XMIN,XMAX,YMIN,YMAX
CHARACTER*(*) TITLE(*)
* Import-Export :
* Export :
* Status :
INTEGER STATUS
* Function declarations :
* Global variables :
INCLUDE 'SPLOT_CMN'
* Local constants :
INTEGER SYMBOL(8)/17,13,16,18,4,7,6,12/
* Local variables :
CHARACTER*40 LABEL
REAL TXTSCALE
REAL X1,X2,Y1,Y2
REAL XW1,XW2,YW1,YW2
REAL BXW1,BXW2,BYW1,BYW2
SAVE BXW1,BXW2,BYW1,BYW2
INTEGER ISYM
*-
* Check inherited global status.
IF ( STATUS .NE. SAI__OK ) RETURN
* Set plot position and text size
X1=XMIN
X2=XMAX
Y1=YMIN
Y2=YMAX
IF(NPZ.EQ.2)THEN
IF(PPLOT)THEN ! Mode DM
Y2=0.5*(YMIN+YMAX)
ELSE ! Mode DR
Y1=YMIN+0.33*(YMAX-YMIN)
END IF
ELSEIF (NPZ.GT.2) THEN ! Mode DMR
X2=0.5*(XMIN+XMAX)
Y1=YMIN+0.33*(YMAX-YMIN)
END IF
TXTSCALE=MIN(0.8,2*(X2-X1))
* Now allow margin for labels & title
X1=X1+3.0*TXTSCALE*0.025 ! 3.0 char. heights
Y1=Y1+3.0*TXTSCALE*0.025 ! 3.0 char. heights
Y2=Y2-2.5*TXTSCALE*0.025 ! 2.5 char. heights
* Viewport position
CALL GCB_SETR('POS_X1',X1,STATUS)
CALL GCB_SETR('POS_X2',X2,STATUS)
CALL GCB_SETR('POS_Y1',Y1,STATUS)
CALL GCB_SETR('POS_Y2',Y2,STATUS)
* numeric label character size
CALL GCB_SETR('AXES_SIZE',TXTSCALE,STATUS)
* axis labels
CALL GCB_SETR('XLABEL_SIZE',TXTSCALE,STATUS)
IF (EAXIS) THEN
LABEL='Energy (keV)'
ELSE
LABEL='Channel no.'
END IF
CALL GCB_SETC('XLABEL_TEXT',LABEL,STATUS)
CALL GCB_SETR('YLABEL_SIZE',TXTSCALE,STATUS)
CALL GCB_SETC('YLABEL_TEXT','count s\u-1\d channel\u-1\d',
: STATUS)
IF (LOGAX) THEN
CALL GCB_SETL('XAXIS_LOG',.TRUE.,STATUS)
CALL GCB_SETL('YAXIS_LOG',.TRUE.,STATUS)
END IF
* Title lines
CALL GCB_SETI('TITLE_N',2,STATUS)
CALL GCB_SET1C('TITLE_TEXT',1,2,TITLE,STATUS)
CALL GCB_SET1R('TITLE_SIZE',1,1,TXTSCALE,STATUS)
CALL GCB_SET1R('TITLE_SIZE',2,1,TXTSCALE,STATUS)
CALL GCB_SET1C('TITLE_JUST',1,1,'L',STATUS)
CALL GCB_SET1C('TITLE_JUST',2,1,'L',STATUS)
* if several datasets to be overlayed on same axes need overall max/min
IF (OVERLAY) THEN
CALL GFX_DEF1DWNDV(NVAL,%val(PXPTR),%val(DDPTR),%val(DVPTR),
: XW1,XW2,YW1,YW2,STATUS)
IF (IDS.EQ.1) THEN
BXW1=XW1
BXW2=XW2
BYW1=YW1
BYW2=YW2
ELSE
BXW1=MIN(XW1,BXW1)
BXW2=MAX(XW2,BXW2)
BYW1=MIN(YW1,BYW1)
BYW2=MAX(YW2,BYW2)
END IF
IF (IDS.EQ.NDS) THEN
CALL GCB_SETR('XAXIS_LO',BXW1,STATUS)
CALL GCB_SETR('XAXIS_HI',BXW2,STATUS)
CALL GCB_SETR('YAXIS_LO',BYW1,STATUS)
CALL GCB_SETR('YAXIS_HI',BYW2,STATUS)
END IF
* set plotting symbol for each dataset
ISYM=MOD(IDS,8)
IF (ISYM.EQ.0) THEN
ISYM=8
END IF
CALL GCB_SETL('POINT_FLAG',.TRUE.,STATUS)
CALL GCB_SETI('POINT_SYMBOL',SYMBOL(ISYM),STATUS)
CALL GCB_SETL('ERR_FLAG',.TRUE.,STATUS)
END IF
END
*+
SUBROUTINE SPLOT_DPLOT_COPY(OFID,IP,IG,IDS,SAVE,STATUS)
* Description :
* Type definitions :
IMPLICIT NONE
* Global constants :
INCLUDE 'SAE_PAR'
INCLUDE 'FIT_PAR'
* Global variables :
INCLUDE 'SPLOT_CMN'
* Structure definitions :
* Import :
INTEGER OFID
INTEGER IP,IG,IDS
LOGICAL SAVE
* Import-Export :
* Export :
* Status :
INTEGER STATUS
* Function declarations :
INTEGER CHR_LEN
* Local constants :
* Local variables :
CHARACTER*80 OVLY
CHARACTER*80 UNITS
INTEGER NDFS(NDSMAX)
INTEGER NDFID,XPTR,WPTR
INTEGER I,L
*-
* Check inherited global status.
IF ( STATUS .NE. SAI__OK ) RETURN
* Axis in channels or energy
IF ( EAXIS ) THEN
XPTR = PXPTR
WPTR = PWPTR
ELSE
XPTR = DXPTR
WPTR = DWPTR
END IF
* Create the data graph
CALL SPLOT_CRENDF( SAVE, 'channel space data',
: OFID, IG, NVAL, DDPTR, DVPTR, DQPTR,
: XPTR, WPTR,
: 'Channel/Data space - data',
: NDFID, STATUS )
* Copy labels
CALL BDI_COPY( DSID, 'Label,Units,Title', NDFID, ' ', STATUS )
* Modify units to take account of normalisation
CALL BDI_GET0C( DSID, 'Units', UNITS, STATUS )
L = CHR_LEN(UNITS) + 1
UNITS(L:) = '/channel'
CALL BDI_PUT0C( NDFID, 'Units', UNITS, STATUS )
* Ancillary info
IF ( SAVE ) THEN
CALL UDI_COPANC( DSID, 'grf', NDFID, STATUS )
END IF
* Release graph
CALL SPLOT_RELNDF( NDFID, STATUS )
* Cancel attributes that not relevant to other data
CALL GCB_CAN( 'TITLE_N,YLABEL_TEXT', STATUS )
CALL GCB_CAN1( 'TITLE_SIZE,TITLE_TEXT,TITLE_JUST', 1, 2, STATUS )
IF ( OVERLAY ) THEN
CALL GCB_CAN( 'POINT_FLAG,POINT_SYMBOL', STATUS )
CALL GCB_CAN( 'ERR_FLAG', STATUS )
CALL GCB_CAN( 'XAXIS_LO,XAXIS_HI,YAXIS_LO,YAXIS_HI', STATUS )
END IF
* Create the fit graph
CALL SPLOT_CRENDF( SAVE, 'channel space fit',
: OFID, IG+1, NVAL, DFDPTR, 0, 0,
: XPTR, WPTR,
: 'Channel/Data space - fit',
: NDFID, STATUS )
* If single point then set marker to '*' for model point
IF ( NVAL .EQ. 1 ) THEN
CALL GCB_SETL('POINT_FLAG',.TRUE.,STATUS)
CALL GCB_SETI('POINT_SYMBOL',3,STATUS)
END IF
* Datasets overlayed - record NDF number
IF (OVERLAY) THEN
NDFS(IDS)=IG
* last dataset - construct overlays
IF (IDS.EQ.NDS) THEN
OVLY=' '
L=0
DO I=1,NDS-1
L=L+6
WRITE(OVLY(L-5:L),'(2(I2,A1))') NDFS(I),',',NDFS(I)+1,','
END DO
L=L+1
WRITE(OVLY(L:L+1),'(I2)') IG+1
CALL GMI_SETPLOT( OFID, IP, IG, OVLY, STATUS )
END IF
ELSE
* set overlay
WRITE(OVLY,'(I2)') IG+1
CALL GMI_SETPLOT( OFID, IP, IG, OVLY, STATUS )
END IF
* Release graph
CALL SPLOT_RELNDF( NDFID, STATUS )
CALL GCB_CAN( 'POINT_FLAG,POINT_SYMBOL', STATUS )
CALL GCB_CAN( 'XAXIS_LOG,YAXIS_LOG,YAXIS_LO,YAXIS_HI', STATUS )
END
*+
SUBROUTINE SPLOT_DRPLOT_SETUP(IDS,NPZ,XMIN,XMAX,YMIN,YMAX,STATUS)
* Description :
* Type definitions :
IMPLICIT NONE
* Global constants :
INCLUDE 'SAE_PAR'
INCLUDE 'FIT_PAR'
* Global variables :
INCLUDE 'SPLOT_CMN'
* Structure definitions :
* Import :
INTEGER IDS,NPZ
REAL XMIN,XMAX,YMIN,YMAX
* Import-Export :
* Export :
* Status :
INTEGER STATUS
* Function declarations :
* Local constants :
INTEGER SYMBOL(8)/17,13,16,18,4,7,6,12/
* Local variables :
CHARACTER*20 LABEL
REAL X1,X2,Y1,Y2
REAL TXTSCALE
REAL XW1,XW2,YW1,YW2
REAL BXW1,BXW2,BYW1,BYW2
SAVE BXW1,BXW2,BYW1,BYW2
INTEGER ISYM
*-
* Check inherited global status.
IF ( STATUS .NE. SAI__OK ) RETURN
* Calculate residuals
CALL SPLOT_RESID(NVAL,%val(DDPTR),%val(DFDPTR),%val(DRDPTR),
: STATUS)
* Set plot position
X1=XMIN
X2=XMAX
Y1=YMIN
Y2=YMAX
IF(NPZ.EQ.2)THEN ! Mode DR
Y2=YMIN+0.33*(YMAX-YMIN)
ELSEIF (NPZ.GT.2) THEN ! Mode DMR
X2=0.5*(XMIN+XMAX)
Y2=YMIN+0.33*(YMAX-YMIN)
END IF
TXTSCALE=MIN(0.8,2*(X2-X1))
* now allow margin for labels & title
X1=X1+3.0*TXTSCALE*0.025 ! 3.0 char. heights
Y1=Y1+3.0*TXTSCALE*0.025 ! 3.0 char. heights
Y2=Y2-2.5*TXTSCALE*0.025 ! 2.5 char. heights
CALL GCB_SETR('POS_X1',X1,STATUS)
CALL GCB_SETR('POS_X2',X2,STATUS)
CALL GCB_SETR('POS_Y1',Y1,STATUS)
CALL GCB_SETR('POS_Y2',Y2,STATUS)
CALL GCB_SETR('AXES_SIZE',TXTSCALE,STATUS)
* axis labels
CALL GCB_SETR('XLABEL_SIZE',TXTSCALE,STATUS)
IF (EAXIS) THEN
LABEL='Energy (keV)'
ELSE
LABEL='Channel no.'
END IF
CALL GCB_SETC('XLABEL_TEXT',LABEL,STATUS)
CALL GCB_SETR('YLABEL_SIZE',TXTSCALE,STATUS)
IF (LOGAX) THEN
CALL GCB_SETL('XAXIS_LOG',.TRUE.,STATUS)
END IF
CALL GCB_SETI('TITLE_N',1,STATUS)
CALL GCB_SET1C('TITLE_TEXT',1,1,'Residuals (data-model)',STATUS)
CALL GCB_SET1R('TITLE_SIZE',1,1,TXTSCALE,STATUS)
* if several datasets to be overlayed on same axes need overall max/min
IF (OVERLAY) THEN
CALL GFX_DEF1DWNDV(NVAL,%val(PXPTR),%val(DRDPTR),%val(DVPTR),
: XW1,XW2,YW1,YW2,STATUS)
IF (IDS.EQ.1) THEN
BXW1=XW1
BXW2=XW2
BYW1=YW1
BYW2=YW2
ELSE
BXW1=MIN(XW1,BXW1)
BXW2=MAX(XW2,BXW2)
BYW1=MIN(YW1,BYW1)
BYW2=MAX(YW2,BYW2)
END IF
IF (IDS.EQ.NDS) THEN
CALL GCB_SETR('XAXIS_LO',BXW1,STATUS)
CALL GCB_SETR('XAXIS_HI',BXW2,STATUS)
CALL GCB_SETR('YAXIS_LO',BYW1,STATUS)
CALL GCB_SETR('YAXIS_HI',BYW2,STATUS)
END IF
* set plotting symbol for each dataset
ISYM=MOD(IDS,8)
IF (ISYM.EQ.0) THEN
ISYM=8
END IF
CALL GCB_SETL('POINT_FLAG',.TRUE.,STATUS)
CALL GCB_SETI('POINT_SYMBOL',SYMBOL(ISYM),STATUS)
CALL GCB_SETL('ERR_FLAG',.TRUE.,STATUS)
END IF
END
*+
SUBROUTINE SPLOT_DRPLOT_COPY(OFID,IP,IG,IDS,SAVE,STATUS)
* Description :
* Type definitions :
IMPLICIT NONE
* Global constants :
INCLUDE 'SAE_PAR'
INCLUDE 'FIT_PAR'
* Global variables :
INCLUDE 'SPLOT_CMN'
* Structure definitions :
* Import :
INTEGER OFID,IP,IG,IDS
LOGICAL SAVE
* Import-Export :
* Export :
* Status :
INTEGER STATUS
* Function declarations :
* Local constants :
* Local variables :
CHARACTER*80 OVLY
INTEGER NDFS(NDSMAX)
INTEGER NDFID,XPTR,WPTR
INTEGER I,L
*-
* Check inherited global status.
IF ( STATUS .NE. SAI__OK ) RETURN
* Axis in channels or energy
IF (EAXIS) THEN
XPTR=PXPTR
WPTR=PWPTR
ELSE
XPTR=DXPTR
WPTR=DWPTR
END IF
* Create the residuals graph
CALL SPLOT_CRENDF( SAVE, 'channel space residuals',
: OFID, IG, NVAL, DRDPTR, DVPTR, DQPTR,
: XPTR, WPTR,
: 'Channel/Data space - residuals',
: NDFID, STATUS )
* Add title
CALL BDI_PUT0C( NDFID, 'Title', 'Residuals (data-model)', STATUS )
* datasets overlayed - record NDF number
IF (OVERLAY) THEN
NDFS(IDS)=IG
* last dataset - construct overlays
IF (IDS.EQ.NDS) THEN
OVLY=' '
L=0
DO I=1,NDS-1
L=L+3
WRITE(OVLY(L-2:L),'(I2,A1)') NDFS(I),','
END DO
OVLY(L:L)=' '
CALL GMI_SETPLOT( OFID, IP, IG, OVLY, STATUS )
END IF
ELSE
CALL GMI_SETPLOT( OFID, IP, IG, ' ', STATUS )
END IF
* Release graph
CALL SPLOT_RELNDF( NDFID, STATUS )
CALL GCB_CAN( 'XAXIS_LOG,YAXIS_LOG', STATUS )
IF ( OVERLAY ) THEN
CALL GCB_CAN( 'POINT_FLAG,POINT_SYMBOL', STATUS )
CALL GCB_CAN( 'ERR_FLAG', STATUS )
CALL GCB_CAN( 'XAXIS_LO,XAXIS_HI,YAXIS_LO,YAXIS_HI', STATUS )
END IF
END
*+
SUBROUTINE SPLOT_PPLOT_SETUP(IDS,NPZ,XMIN,XMAX,YMIN,YMAX,
: TITLE,STATUS)
* Description :
* Type definitions :
IMPLICIT NONE
* Global constants :
INCLUDE 'SAE_PAR'
INCLUDE 'FIT_PAR'
* Global variables :
INCLUDE 'SPLOT_CMN'
* Structure definitions :
* Import :
INTEGER IDS,NPZ
REAL XMIN,XMAX,YMIN,YMAX
CHARACTER*(*) TITLE(*)
* Import-Export :
* Export :
* Status :
INTEGER STATUS
* Function declarations :
* Local constants :
INTEGER SYMBOL(8)/17,13,16,18,4,7,6,12/
* Local variables :
REAL X1,X2,Y1,Y2
REAL TXTSCALE
REAL XW1,XW2,YW1,YW2
REAL BXW1,BXW2,BYW1,BYW2
SAVE BXW1,BXW2,BYW1,BYW2
INTEGER ISYM
LOGICAL ERR
*-
* Check inherited global status.
IF ( STATUS .NE. SAI__OK ) RETURN
ERR = .TRUE.
* Set up log:log flag
CALL GCB_SETL('XAXIS_LOG',.TRUE.,STATUS)
CALL GCB_SETL('YAXIS_LOG',.TRUE.,STATUS)
* Set plot position
X1=XMIN ! Mode P
X2=XMAX
Y1=YMIN
Y2=YMAX
IF(NPZ.EQ.2)THEN
IF(DPLOT)THEN ! Mode CP
Y1=0.5*(YMIN+YMAX)
ELSE ! Mode PR
Y1=YMIN+0.33*(YMAX-YMIN)
END IF
ELSEIF (NPZ.GT.2) THEN ! Mode CPR
X1=0.5*(XMIN+XMAX)
Y1=YMIN+0.33*(YMAX-YMIN)
END IF
TXTSCALE=MIN(0.8,2*(X2-X1))
* now allow margin for labels & title
X1=X1+3.0*TXTSCALE*0.025 ! 3.0 char. heights
Y1=Y1+3.0*TXTSCALE*0.025 ! 3.0 char. heights
Y2=Y2-2.5*TXTSCALE*0.025 ! 2.5 char. heights
CALL GCB_SETR('POS_X1',X1,STATUS)
CALL GCB_SETR('POS_X2',X2,STATUS)
CALL GCB_SETR('POS_Y1',Y1,STATUS)
CALL GCB_SETR('POS_Y2',Y2,STATUS)
CALL GCB_SETR('AXES_SIZE',TXTSCALE,STATUS)
CALL GCB_SETC('XLABEL_TEXT','Energy (keV)',STATUS)
CALL GCB_SETR('XLABEL_SIZE',TXTSCALE,STATUS)
CALL GCB_SETC('YLABEL_TEXT',
: 'photon cm\u-2\d s\u-1\d keV\u-1\d',STATUS)
CALL GCB_SETR('YLABEL_SIZE',TXTSCALE,STATUS)
CALL GCB_SET1C('TITLE_TEXT',1,2,TITLE,STATUS)
CALL GCB_SETI('TITLE_N',2,STATUS)
CALL GCB_SET1R('TITLE_SIZE',1,1,TXTSCALE,STATUS)
CALL GCB_SET1R('TITLE_SIZE',2,1,TXTSCALE,STATUS)
CALL GCB_SET1C('TITLE_JUST',1,1,'L',STATUS)
CALL GCB_SET1C('TITLE_JUST',2,1,'L',STATUS)
* if several datasets to be overlayed on same axes need overall max/min
IF (OVERLAY) THEN
CALL GFX_DEF1DWNDV(NVAL,%val(PXPTR),%val(PDPTR),%val(PVPTR),
: XW1,XW2,YW1,YW2,STATUS)
IF (IDS.EQ.1) THEN
BXW1=XW1
BXW2=XW2
BYW1=YW1
BYW2=YW2
ELSE
BXW1=MIN(XW1,BXW1)
BXW2=MAX(XW2,BXW2)
BYW1=MIN(YW1,BYW1)
BYW2=MAX(YW2,BYW2)
END IF
IF (IDS.EQ.NDS) THEN
CALL GCB_SETR('XAXIS_LO',BXW1,STATUS)
CALL GCB_SETR('XAXIS_HI',BXW2,STATUS)
CALL GCB_SETR('YAXIS_LO',BYW1,STATUS)
CALL GCB_SETR('YAXIS_HI',BYW2,STATUS)
END IF
* set plotting symbol for each dataset
ISYM=MOD(IDS,8)
IF (ISYM.EQ.0) THEN
ISYM=8
END IF
CALL GCB_SETL('POINT_FLAG',.TRUE.,STATUS)
CALL GCB_SETI('POINT_SYMBOL',SYMBOL(ISYM),STATUS)
CALL GCB_SETL('ERR_FLAG',.TRUE.,STATUS)
END IF
END
*+
SUBROUTINE SPLOT_PPLOT_COPY(OFID,IP,IG,IDS,SAVE,STATUS)
* Description :
* Type definitions :
IMPLICIT NONE
* Global constants :
INCLUDE 'SAE_PAR'
INCLUDE 'FIT_PAR'
* Global variables :
INCLUDE 'SPLOT_CMN'
* Structure definitions :
* Import :
INTEGER OFID,IP,IG,IDS
LOGICAL SAVE
* Import-Export :
* Export :
* Status :
INTEGER STATUS
* Function declarations :
* Local constants :
* Local variables :
CHARACTER*80 OVLY
INTEGER NDFS(NDSMAX)
INTEGER NDFID
INTEGER I,L
*-
* Check inherited global status.
IF ( STATUS .NE. SAI__OK ) RETURN
* Model space data (i.e. transformed back using fitted model)
CALL SPLOT_CRENDF( SAVE, 'photon space data',
: OFID, IG, NVAL, PDPTR, PVPTR, DQPTR,
: PXPTR, PWPTR,
: 'Photon space - data',
: NDFID, STATUS )
* Copy data label and create new units
CALL BDI_COPY( DSID, 'Label', NDFID, ' ', STATUS )
CALL BDI_PUT0C( NDFID, 'Units', 'photon/(cm**2*sec*keV)', STATUS )
* Create axis label and units
CALL BDI_AXPUT0C( NDFID, 1, 'Label', 'Energy', STATUS )
CALL BDI_AXPUT0C( NDFID, 1, 'Units', 'keV', STATUS )
* Release graph
CALL SPLOT_RELNDF( NDFID, STATUS )
CALL GCB_CAN( 'TITLE_N', STATUS )
CALL GCB_CAN1( 'TITLE_SIZE,TITLE_TEXT,TITLE_JUST', 1, 2, STATUS )
CALL GCB_CAN( 'YLABEL_TEXT', STATUS )
CALL GCB_CAN( 'YAXIS_LO,YAXIS_HI', STATUS )
IF ( OVERLAY ) THEN
CALL GCB_CAN( 'POINT_FLAG,POINT_SYMBOL', STATUS )
CALL GCB_CAN( 'ERR_FLAG', STATUS )
CALL GCB_CAN( 'XAXIS_LO,XAXIS_HI', STATUS )
END IF
* Model space fit (for overlay)
CALL SPLOT_CRENDF( SAVE, 'photon space fit',
: OFID, IG+1, NFIT, PFDPTR, 0, 0,
: PFXPTR, 0,
: 'Photon space - fit',
: NDFID, STATUS )
* Energy axis
CALL BDI_AXPUT0C( NDFID, 1, 'Label', 'Energy', STATUS )
CALL BDI_AXPUT0C( NDFID, 1, 'Units', 'keV', STATUS )
* Release graph
CALL SPLOT_RELNDF( NDFID, STATUS )
* Datasets overlaid - record NDF number
IF ( OVERLAY ) THEN
NDFS(IDS)=IG
* last dataset - construct overlays
IF (IDS.EQ.NDS) THEN
OVLY=' '
L=0
DO I=1,NDS-1
L=L+6
WRITE(OVLY(L-5:L),'(2(I2,A1))') NDFS(I),',',NDFS(I)+1,','
END DO
L=L+1
WRITE(OVLY(L:L+1),'(I2)') IG+1
CALL GMI_SETPLOT( OFID, IP, IG, OVLY, STATUS )
END IF
ELSE
* set overlay
WRITE(OVLY,'(I2)') IG+1
CALL GMI_SETPLOT( OFID, IP, IG, OVLY, STATUS )
END IF
CALL GCB_CAN( 'XAXIS_LOG,YAXIS_LOG', STATUS )
END
*+
SUBROUTINE SPLOT_PRPLOT_SETUP(IDS,NPZ,XMIN,XMAX,YMIN,YMAX,STATUS)
* Description :
* Type definitions :
IMPLICIT NONE
* Global constants :
INCLUDE 'SAE_PAR'
INCLUDE 'QUAL_PAR'
INCLUDE 'FIT_PAR'
* Global variables :
INCLUDE 'SPLOT_CMN'
* Structure definitions :
* Import :
INTEGER IDS,NPZ
REAL XMIN,XMAX,YMIN,YMAX
* Import-Export :
* Export :
* Status :
INTEGER STATUS
* Function declarations :
* Local constants :
INTEGER SYMBOL(8)/17,13,16,18,4,7,6,12/
* Local variables :
REAL X1,X2,Y1,Y2
REAL TXTSCALE
REAL XW1,XW2,YW1,YW2
REAL BXW1,BXW2,BYW1,BYW2
SAVE BXW1,BXW2,BYW1,BYW2
INTEGER ISYM
*-
* Check inherited global status.
IF ( STATUS .NE. SAI__OK ) RETURN
IF ( QOK ) THEN
CALL ARR_COP1B( NVAL, %VAL(DQPTR), %VAL(PRQPTR), STATUS )
ELSE
CALL ARR_INIT1B( QUAL__GOOD, NVAL, %VAL(PRQPTR), STATUS )
END IF
* Subtract fit from data values & divide by sqrt(variance)
CALL SPLOT_MODRES( NVAL, %VAL(PDPTR), %VAL(PIDPTR),
: %VAL(PVPTR), %VAL(PRDPTR), %VAL(PRVPTR),
: %VAL(PRQPTR), STATUS )
* Set plot position
IF ( NPZ .EQ. 2 ) THEN ! Mode MR
X1 = XMIN
X2 = XMAX
Y1 = YMIN
Y2 = YMIN+0.33*(YMAX-YMIN)
ELSE ! Mode DMR
X1 = 0.5*(XMIN+XMAX)
X2 = XMAX
Y1 = YMIN
Y2 = YMIN+0.33*(YMAX-YMIN)
END IF
TXTSCALE = MIN(0.8,2*(X2-X1))
* Now allow margin for labels & title
X1=X1+3.0*TXTSCALE*0.025 ! 3.0 char. heights
Y1=Y1+3.0*TXTSCALE*0.025 ! 3.0 char. heights
Y2=Y2-2.5*TXTSCALE*0.025 ! 2.5 char. heights
CALL GCB_SETR('POS_X1',X1,STATUS)
CALL GCB_SETR('POS_X2',X2,STATUS)
CALL GCB_SETR('POS_Y1',Y1,STATUS)
CALL GCB_SETR('POS_Y2',Y2,STATUS)
CALL GCB_SETR('AXES_SIZE',TXTSCALE,STATUS)
CALL GCB_SETR('XLABEL_SIZE',TXTSCALE,STATUS)
CALL GCB_SETC('XLABEL_TEXT','Energy (keV)',STATUS)
CALL GCB_SETR('YLABEL_SIZE',TXTSCALE,STATUS)
CALL GCB_SETC('YLABEL_TEXT','Significance',STATUS)
CALL GCB_SET1C('TITLE_TEXT',1,1,
: 'Normalised residuals (data-model)/sigma',STATUS)
CALL GCB_SET1R('TITLE_SIZE',1,1,TXTSCALE,STATUS)
CALL GCB_SETI('TITLE_N',1,STATUS)
CALL GCB_SETL('XAXIS_LOG',.TRUE.,STATUS)
CALL GCB_SETL('YAXIS_LOG',.FALSE.,STATUS)
* if several datasets to be overlayed on same axes need overall max/min
IF (OVERLAY) THEN
CALL GFX_DEF1DWNDV(NVAL,%val(PXPTR),%val(PRDPTR),%val(PRVPTR),
: XW1,XW2,YW1,YW2,STATUS)
IF (IDS.EQ.1) THEN
BXW1=XW1
BXW2=XW2
BYW1=YW1
BYW2=YW2
ELSE
BXW1=MIN(XW1,BXW1)
BXW2=MAX(XW2,BXW2)
BYW1=MIN(YW1,BYW1)
BYW2=MAX(YW2,BYW2)
END IF
IF (IDS.EQ.NDS) THEN
CALL GCB_SETR('XAXIS_LO',BXW1,STATUS)
CALL GCB_SETR('XAXIS_HI',BXW2,STATUS)
CALL GCB_SETR('YAXIS_LO',BYW1,STATUS)
CALL GCB_SETR('YAXIS_HI',BYW2,STATUS)
END IF
* Set plotting symbol for each dataset
ISYM=MOD(IDS,8)
IF (ISYM.EQ.0) THEN
ISYM=8
END IF
CALL GCB_SETL('POINT_FLAG',.TRUE.,STATUS)
CALL GCB_SETI('POINT_SYMBOL',SYMBOL(ISYM),STATUS)
CALL GCB_SETL('ERR_FLAG',.TRUE.,STATUS)
END IF
END
*+
SUBROUTINE SPLOT_PRPLOT_COPY(OFID,IP,IG,IDS,SAVE,STATUS)
* Description :
* Type definitions :
IMPLICIT NONE
* Global constants :
INCLUDE 'SAE_PAR'
INCLUDE 'FIT_PAR'
* Global variables :
INCLUDE 'SPLOT_CMN'
* Structure definitions :
* Import :
INTEGER OFID,IP,IG,IDS
LOGICAL SAVE
* Import-Export :
* Export :
* Status :
INTEGER STATUS
* Function declarations :
* Local constants :
* Local variables :
CHARACTER*80 OVLY
INTEGER NDFS(NDSMAX)
INTEGER NDFID
INTEGER I,L
*-
* Check inherited global status.
IF ( STATUS .NE. SAI__OK ) RETURN
* Create residuals graph
CALL SPLOT_CRENDF( SAVE, 'photon space residuals',
: OFID, IG, NVAL, PRDPTR, PRVPTR, PRQPTR,
: PXPTR, PWPTR,
: 'Photon space - residuals',
: NDFID, STATUS )
CALL BDI_PUT0C( NDFID, 'Title',
: 'Normalised residuals (data-model)/sigma', STATUS )
CALL BDI_PUT0C( NDFID, 'Label', 'Significance', STATUS )
CALL BDI_AXPUT0C( NDFID, 1, 'Label', 'Energy', STATUS )
CALL BDI_AXPUT0C( NDFID, 1, 'Units', 'keV', STATUS )
* datasets overlayed - record NDF number
IF (OVERLAY) THEN
NDFS(IDS)=IG
* last dataset - construct overlays
IF (IDS.EQ.NDS) THEN
OVLY=' '
L=0
DO I=1,NDS-1
L=L+3
WRITE(OVLY(L-2:L),'(I2,A1)') NDFS(I),','
END DO
OVLY(L:L)=' '
CALL GMI_SETPLOT( OFID, IP, IG, OVLY, STATUS )
END IF
ELSE
CALL GMI_SETPLOT( OFID, IP, IG, ' ', STATUS )
END IF
* Release graph
CALL SPLOT_RELNDF( NDFID, STATUS )
CALL GCB_CAN( 'TITLE_N', STATUS )
CALL GCB_CAN1( 'TITLE_SIZE,TITLE_TEXT', 1, 2, STATUS )
CALL GCB_CAN( 'XAXIS_LOG,YAXIS_LOG', STATUS )
IF ( OVERLAY ) THEN
CALL GCB_CAN( 'POINT_FLAG,POINT_SYMBOL', STATUS )
CALL GCB_CAN( 'ERR_FLAG', STATUS )
CALL GCB_CAN( 'XAXIS_LO,XAXIS_HI,YAXIS_LO,YAXIS_HI', STATUS )
END IF
END
*+
SUBROUTINE SPLOT_MPLOT_SETUP(NCHAN,NTERM,MODSPEC,COMPSPEC,STATUS)
* Description :
* Type definitions :
IMPLICIT NONE
* Global constants :
INCLUDE 'SAE_PAR'
INCLUDE 'FIT_PAR'
* Global variables :
* Structure definitions :
* Import :
INTEGER NCHAN,NTERM
CHARACTER*(*) MODSPEC,COMPSPEC(*)
* Import-Export :
* Export :
* Status :
INTEGER STATUS
* Function declarations :
* Global variables :
INCLUDE 'SPLOT_CMN'
* Local constants :
* Local variables :
CHARACTER*80 TEXT
REAL TXTSCALE
REAL X1,X2,Y1,Y2
REAL XW1,XW2,YW1,YW2
REAL XRAN,YRAN
REAL XL,YL,XT,YT,XTL,YTL
REAL XDISP,YSTEP
INTEGER I
INTEGER STYLE
*-
* Check inherited global status.
IF ( STATUS .NE. SAI__OK ) RETURN
* Set plot position and text size
X1 = 0.0
X2 = 0.99
Y1 = 0.0
Y2 = 0.99
TXTSCALE=MIN(0.8,2*(X2-X1))
* Now allow margin for labels & title
X1=X1+3.0*TXTSCALE*0.025 ! 3.0 char. heights
Y1=Y1+3.0*TXTSCALE*0.025 ! 3.0 char. heights
Y2=Y2-2.5*TXTSCALE*0.025 ! 2.5 char. heights
* Viewport position
CALL GCB_SETR( 'POS_X1', X1, STATUS )
CALL GCB_SETR( 'POS_X2', X2, STATUS )
CALL GCB_SETR( 'POS_Y1', Y1, STATUS )
CALL GCB_SETR( 'POS_Y2', Y2, STATUS )
IF ( LOGAX ) THEN
CALL GCB_SETL('XAXIS_LOG',.TRUE.,STATUS)
CALL GCB_SETL('YAXIS_LOG',.TRUE.,STATUS)
END IF
* Get axis limits for key
CALL GFX_DEF1DWND(NCHAN,%val(MCXPTR),%val(MTDPTR),
: XW1,XW2,YW1,YW2,STATUS)
* Get positions for key
XRAN=XW2-XW1
YRAN=YW2-YW1
XDISP=0.1*XRAN
XL=XW1+0.8*XRAN
YL=YW1+0.9*YRAN
XT=XW1+0.91*XRAN
YT=YL-0.5*TXTSCALE*0.025*YRAN
YSTEP=0.05*YRAN
c CALL GCB_SETI('SHAPE_N',NTERM+1,STATUS)
c CALL GCB_SET1C('SHAPE_TYPE',1,1,'VECTOR',STATUS)
c CALL GCB_SET1R('SHAPE_DATA1',1,1,XDISP,STATUS)
c CALL GCB_SET1R('SHAPE_DATA2',1,1,0.0,STATUS)
c CALL GCB_SET1R('SHAPE_X',1,1,XL,STATUS)
c CALL GCB_SET1R('SHAPE_Y',1,1,YL,STATUS)
c CALL GCB_SETI('NOTE_N',NTERM+1,STATUS)
IF (LOGAX) THEN
XTL=10.0**XT
YTL=10.0**YT
c CALL GCB_SET1R('NOTE_X',1,1,XTL,STATUS)
c CALL GCB_SET1R('NOTE_Y',1,1,YTL,STATUS)
ELSE
c CALL GCB_SET1R('NOTE_X',1,1,XT,STATUS)
c CALL GCB_SET1R('NOTE_Y',1,1,YT,STATUS)
END IF
c CALL GCB_SET1C('NOTE_TEXT',1,1,'Total',STATUS)
c CALL GCB_SET1R('NOTE_SIZE',1,1,TXTSCALE,STATUS)
STYLE=0
DO I=1,NTERM
STYLE=STYLE+1
STYLE=STYLE/5 + MOD(STYLE,5) + 1
c CALL GCB_SET1C('SHAPE_TYPE',I+1,1,'VECTOR',STATUS)
c CALL GCB_SET1R('SHAPE_DATA1',I+1,1,XDISP,STATUS)
c CALL GCB_SET1R('SHAPE_DATA2',I+1,1,0.0,STATUS)
c CALL GCB_SET1I('SHAPE_STYLE',I+1,1,STYLE,STATUS)
c CALL GCB_SET1R('SHAPE_X',I+1,1,XL,STATUS)
YL=YL-YSTEP
c CALL GCB_SET1R('SHAPE_Y',I+1,1,YL,STATUS)
c CALL GCB_SET1C('NOTE_TEXT',I+1,1,COMPSPEC(I),STATUS)
c CALL GCB_SET1R('NOTE_SIZE',I+1,1,TXTSCALE,STATUS)
YT=YT-YSTEP
IF (LOGAX) THEN
YTL=10.0**YT
c CALL GCB_SET1R('NOTE_X',I+1,1,XTL,STATUS)
c CALL GCB_SET1R('NOTE_Y',I+1,1,YTL,STATUS)
ELSE
c CALL GCB_SET1R('NOTE_X',I+1,1,XT,STATUS)
c CALL GCB_SET1R('NOTE_Y',I+1,1,YT,STATUS)
END IF
END DO
* Numeric label character size
CALL GCB_SETR('AXES_SIZE',TXTSCALE,STATUS)
* Axis labels
CALL GCB_SETR('XLABEL_SIZE',TXTSCALE,STATUS)
TEXT='Energy (keV)'
CALL GCB_SETC('XLABEL_TEXT',TEXT,STATUS)
CALL GCB_SETR('YLABEL_SIZE',TXTSCALE,STATUS)
TEXT='photon s\u-1\d cm\u-2\d keV\u-1\d'
CALL GCB_SETC('YLABEL_TEXT',TEXT,STATUS)
* Title lines
CALL GCB_SETI('TITLE_N',1,STATUS)
TEXT = 'Additive components for model: '//MODSPEC
CALL GCB_SET1C('TITLE_TEXT',1,1,TEXT,STATUS)
CALL GCB_SET1R('TITLE_SIZE',1,1,TXTSCALE,STATUS)
CALL GCB_SET1C('TITLE_JUST',1,1,'L',STATUS)
END
*+
SUBROUTINE SPLOT_MPLOT_COPY(OFID,NCHAN,NTERM,COMPSPEC,SAVE,STATUS)
* Description :
* Type definitions :
IMPLICIT NONE
* Global constants :
INCLUDE 'SAE_PAR'
INCLUDE 'FIT_PAR'
* Global variables :
INCLUDE 'SPLOT_CMN'
* Structure definitions :
* Import :
INTEGER OFID,NCHAN,NTERM
CHARACTER*(*) COMPSPEC(*)
LOGICAL SAVE
* Import-Export :
* Export :
* Status :
INTEGER STATUS
* Function declarations :
* Local constants :
* Local variables :
CHARACTER*80 OVLY
INTEGER NDFID
INTEGER IG
INTEGER I
INTEGER STYLE
INTEGER L
*-
* Check inherited global status.
IF ( STATUS .NE. SAI__OK ) RETURN
* locate correct NDF and write index entry
IG = 1
CALL SPLOT_CRENDF( SAVE, 'total model',
: OFID, IG, NCHAN, MTDPTR, 0, 0,
: MCXPTR, 0, 'Total model',
: NDFID, STATUS )
* Release graph
CALL SPLOT_RELNDF( NDFID, STATUS )
* Now each component model
L=1
DO I=1,NTERM
IG=IG+1
CALL SPLOT_CRENDF( SAVE, 'component '//COMPSPEC(I),
: OFID, IG, NCHAN, MCDPTR(I), 0, 0,
: MCXPTR, 0,
: COMPSPEC(I), NDFID, STATUS )
STYLE = STYLE+1
STYLE = STYLE/5 + MOD(STYLE,5) + 1
CALL GCB_SETI('POLY_STYLE',STYLE,STATUS)
* Release graph
CALL SPLOT_RELNDF( NDFID, STATUS )
WRITE(OVLY(L:),'(I2)') IG
IF (I.LT.NTERM) THEN
L=L+2
OVLY(L:L)=','
L=L+1
END IF
END DO
* Cancel attributes that not relevant to other data
CALL GCB_CAN( 'SHAPE_N,NOTE_N', STATUS )
CALL GCB_CAN1( 'SHAPE_TYPE,SHAPE_X,SHAPE_Y,'/
: /'SHAPE_DATA1,SHAPE_DATA2,SHAPE_STYLE',1,NTERM+1,STATUS)
CALL GCB_CAN1( 'NOTE_TEXT,NOTE_X,NOTE_Y', 1, NTERM+1, STATUS )
CALL GCB_CAN( 'POLY_STYLE', STATUS )
CALL GCB_CAN( 'TITLE_N', STATUS )
CALL GCB_CAN1( 'TITLE_SIZE,TITLE_TEXT,TITLE_JUST', 1, 1, STATUS )
CALL GCB_CAN( 'YLABEL_TEXT', STATUS )
CALL GCB_CAN( 'XAXIS_LOG,YAXIS_LOG',STATUS)
CALL GCB_CAN( 'YAXIS_LO,YAXIS_HI', STATUS )
* Set overlay
CALL GMI_SETPLOT( OFID, 1, 1, OVLY, STATUS )
END
*+
SUBROUTINE SPLOT_MPLOT_ADDCOMP(NCHAN,COMP,TOT,STATUS)
* Description :
* Type definitions :
IMPLICIT NONE
* Global constants :
INCLUDE 'SAE_PAR'
* Global variables :
* Structure definitions :
* Import :
INTEGER NCHAN
REAL COMP(*)
* Import-Export :
* Export :
REAL TOT(*)
* Status :
INTEGER STATUS
* Function declarations :
* Local constants :
* Local variables :
INTEGER I
*-
IF (STATUS.EQ.SAI__OK) THEN
DO I=1,NCHAN
TOT(I)=TOT(I)+COMP(I)
END DO
END IF
END
*+
SUBROUTINE SPLOT_MPLOT_EAXIS(NCHAN,EL,EU,ELB,EUB,X,STATUS)
* Description :
* Type definitions :
IMPLICIT NONE
* Global constants :
INCLUDE 'SAE_PAR'
* Import :
INTEGER NCHAN
REAL EL,EU
* Import-Export :
* Export :
REAL ELB(*),EUB(*)
REAL X(*)
* Status :
INTEGER STATUS
* Function declarations :
* Local constants :
* Local variables :
REAL WID,HWID
INTEGER ICHAN
*-
IF (STATUS.EQ.SAI__OK) THEN
WID=(EU-EL)/REAL(NCHAN)
HWID=WID/2.0
ELB(1)=EL
EUB(1)=EL+WID
X(1)=EL+HWID
DO ICHAN=2,NCHAN
ELB(ICHAN)=ELB(ICHAN-1)+WID
EUB(ICHAN)=ELB(ICHAN)+WID
X(ICHAN)=ELB(ICHAN)+HWID
END DO
END IF
END
*+
SUBROUTINE SPLOT_CRENDF( SAVE, ANCE,
: MGID, N, NDAT, DPTR, VPTR, QPTR, XPTR,
: WPTR, NAME, GFID, STATUS )
* Description :
* Type definitions :
IMPLICIT NONE
* Global constants :
INCLUDE 'SAE_PAR'
INCLUDE 'FIT_PAR'
* Import :
LOGICAL SAVE
INTEGER MGID, N, NDAT, DPTR, VPTR, QPTR, XPTR, WPTR
CHARACTER*(*) ANCE,NAME
* Export:
INTEGER GFID
* Global Variables:
INCLUDE 'SPLOT_CMN'
* Status :
INTEGER STATUS
* Local variables :
*-
* Check inherited global status.
IF ( STATUS .NE. SAI__OK ) RETURN
* Annouce this graph
IF ( SAVE ) THEN
CALL MSG_SETI( 'NDF', N )
CALL MSG_PRNT(' Writing '//ANCE//' to NDF ^NDF ...')
END IF
* Locate the graph file object
CALL GMI_LOCNDF( MGID, N, '*', GFID, STATUS )
* Create binned data object and link it to the graph
CALL BDI_LINK( 'Spectrum', 1, NDAT, 'REAL', GFID, STATUS )
* Write the index entry
CALL GMI_PUTINDEX( MGID, N, NAME, STATUS )
* Graph data
CALL BDI_PUT1R( GFID, 'Data', NDAT, %VAL(DPTR), STATUS )
* Graph variance
IF ( VPTR .NE. 0 ) THEN
CALL BDI_PUT1R( GFID, 'Variance', NDAT, %VAL(VPTR), STATUS )
END IF
* Graph quality
IF ( QPTR .NE. 0 ) THEN
CALL BDI_PUT1UB( GFID, 'Quality', NDAT, %VAL(QPTR), STATUS )
CALL BDI_PUT0UB( GFID, 'QualityMask', MASK, STATUS )
END IF
* Axis info
CALL BDI_AXPUT1R( GFID, 1, 'Data', NDAT, %VAL(XPTR), STATUS )
IF ( WPTR .NE. 0 ) THEN
CALL BDI_AXPUT1R( GFID, 1, 'Width', NDAT, %VAL(WPTR), STATUS )
END IF
END
*+
SUBROUTINE SPLOT_RELNDF( NDFID, STATUS )
* Description :
* Type definitions :
IMPLICIT NONE
* Global constants :
INCLUDE 'SAE_PAR'
* Import :
INTEGER NDFID
* Status :
INTEGER STATUS
* Local variables :
*-
* Check inherited global status.
IF ( STATUS .NE. SAI__OK ) RETURN
* Store graphics control
CALL GCB_FSAVE( NDFID, STATUS )
* Release the graph
CALL ADI_ERASE( NDFID, STATUS )
END