subroutine tau_fits C+------------------------------------------------------------------------------ C C ---------------- C T A U _ F I T S C ---------------- C C Description C ----------- C Do fits to lines in a corrected TAURUS cube. C C C Scope of program C ---------------- C - Handles 3-D images only. C - Data array converted to FLOAT. C - Magic values supported. C - Batch execution supported. C - Quality and error arrays not supported. C C C Environment C ----------- C FIGARO C C C Parameters (read or written) C ---------------------------- C NAME_OBS Name of the data cube for analysis. Must be ZXY sorted. C (character)(prompted for). C C ALAM0 Wavelength zeropoint in angstroms. The program first tries C to read this from the cube .MORE.TAURUS structure. C (real)(prompted for) C C DLAM Wavelength increment in angstroms/pixel. The program first tries C to read this from the cube .MORE.TAURUS structure. C (real)(prompted for) C C FSR Free spectral range in angstroms. The program first tries C to read this from the cube .MORE.TAURUS structure. C (real) (prompted for). C C NAME_FIT Name of the output file. (character)(prompted for). C C OPT Menu option. (character)(prompted for) C C BOXSIZE Size of smoothing box in the spatial direction in pixels C (integer)(prompted for) C C NTRIV Number of trivial changes to chi-squared for the fit to be C considered to have converged. (integer)(prompted for) C C DCHI 'Definition' of a trivial fractional change to chi-squared. C (real)(prompted for) C C NITER The maximum number of iterations for any one spectrum. C (integer)(prompted for) C C HOW_GUESS The choice of guessing algorithm. (integer)(prompted for) C C WT_SCHEME The choice of weighting for the data. (integer)(prompted for). C C TYPE_FIT The choice of fitting function. (integer)(prompted for). C C MASK_TRHESH The two extremes in the mask image, between which the spectra C are to be chosen (real[2])(prompted for) C C RNG_AMP, Ranges in amplitude,position, width and background parameter to keep C RNG_POS, from former fits. Used in delimiting which spectra to act on. C RNG_SIG, (real[2])(prompted for) C RNG_BACK C C NAME_MASK Name of file with mask image. (character)(prompted for) C C XTRACT_FIELD Name of file to which to write out an extracted field C (character)(prompted for) C C XTRACT_PAR Number of the parameter which will be written to a separate C file. (integer)(prompted for) C C LAM_CEN Central wavelength for calculation of velocity field. C (real)(prompted for) C C START Starting pixel coords when delimiting a region. C (real[2])(prompted for) C C END Ending pixel coords when delimiting a region. C (real[2])(prompted for) C C REJ_DOM Number of the domain in which one wants to reject fits. C (integer[2])(prompted for) C C NSWEEP Number of sweeps through the frame in NEIGHBOUR in an attempt C to get all of the spectra. (integer)(prompted for) C C BATCH_FIT The type of fitting desired in batch mode. (1) fit, with some C guessing algorithm, (2) nearest neighbour fitting. (integer) C (prompted for) C C TOL Tolerance in percentage of a free spectral range used for cleaning C velocity field. (real)(prompted for) C C VEL_FIELD Name of input/output file for velocity field info. (character) C (prompted for) C C Keywords C -------- C USE_MASK Instruction to use mask image to tell which spectra are to be C fit (prompted for) C C NEW_FILE TRUE if the output file is new. (prompted for) C C C Propagation of data structure C ----------------------------- C Not relevant. C C C Method C ------ C - The IMAGE structure is tested to see if it is 3-dimensional and C is mapped as FLOAT. C - The OUTPUT structure is created to have the same x,y dimensions as C the input cube. The third dimension is the fit parameters (e.g. C peak,centre,width,background etc.) This also has an error C array. The monitor array is pasted into the three extra frames C on top. C - The program does a rough guess for each line depending on which C guessing algorithm the user has chosen. C - The main fitting routine then does the final fit. The fits are C immediately stored away. C - The user can then look at the fits. C C C External functions & subroutines called C --------------------------------------- C Library DSA: C DSA_CLOSE C DSA_CLOSE_STRUCTURE C DSA_DATA_SIZE C DSA_FMTCON C DSA_FREE_LU C DSA_FREE_WORKSPACE C DSA_GET_FLAG_VALUE C DSA_GET_WORKSPACE C DSA_GET_WORK_ARRAY C DSA_INPUT C DSA_MAP_AXIS_DATA C DSA_MAP_DATA C DSA_MAP_ERRORS C DSA_NAMED_INPUT C DSA_NAMED_OUTPUT C DSA_OPEN C DSA_OPEN_TEXT_FILE C DSA_OUTPUT C DSA_POST_PROCESS_FLAGGED_VALUES C DSA_RESHAPE_AXIS C DSA_SET_FLAGGED_VALUES C DSA_SIMPLE_OUTPUT C DSA_TYPESIZE C DSA_UNMAP C DSA_USE_FLAGGED_VALUES C DSA_WRFLUSH C DSA_WRUSER C C Library DYN: C DYN_ELEMENT C C Library GEN: C GEN_CFILL C GEN_FILL C GEN_MOVE C C Library ICH: C ICH_FOLD C ICH_LEN C C Library NDP: C NDP_GET_IMAGE_INFO C NDP_IMAGE_CURSOR C NDP_IMAGE_INDEX C NDP_IMAGE_LUT C NDP_IMAGE_PLOT C NDP_IMAGE_VIEWPORT C NDP_PAR_RDARY C NDP_PGBIN C C Library PAR: C PAR_BATCH C PAR_GIVEN C PAR_RDARY C PAR_RDCHAR C PAR_RDKEY C PAR_RDVAL C PAR_QNUM C PAR_QUEST C PAR_QSTR C C PGPLOT: C PGADVANCE C PGASK C PGBEGIN C PGBOX C PGCURSE C PGDRAW C PGEND C PGLINE C PGMOVE C PGMTEXT C PGPOINT C PGQINF C PGQVP C PGSCH C PGSCI C PGSCR C PGVSIZE C PGVSTAND C PGWINDOW C C Library TAU: C MRQMIN (TAU_PRESS) C TAU_CAL_DTA_NAME C TAU_CAL_RDNUM C TAU_CONVERGE C C Library VAR: C VAR_GETCHR C VAR_GETNUM C VAR_SETCHR C VAR_SETNUM C C C Internal subroutines called C --------------------------- C ACCEPT_FIT C BATCH_OPT C CAUCHY C CAUCHY_EVAL C CLEANVEL C CLEAR_SCREEN C CONV_PARS C CURSE_LIMITS C DISP C EDIT C FAKE_LIM C FIDDLE C FILLIN_FIELD C FIND_CELL C FIND_DOMAIN C FIND_GUESS C FITT C FIT_SPEC C GAUSS C GAUSS_EVAL C GUESS C GUESS_BY_HAND C GUESS_MOMS C GUESS_OTHER C GUESS_PREVIOUS C INFO C LIMITS C MINMAX C NEIGHBOUR C OPTIONS C OUTPUT C PLOT_PANEL C PREPARE_DATA C PREPARE_MON C PRE_PLOT C RD_MON C REJECT C SET_MON C SUB_FIT C TCURSOR C TWO_D_PLOT C TYPE_FIT C VELOCITY_FIELD C WEIGHTS C WHICH_SPECTRUM C WIPE_OUT C ZAP_FIT C C C INCLUDE statements C ------------------ C INCLUDE 'DYNAMIC_MEMORY' C INCLUDE 'NDP_NUMERIC_RANGES' C INCLUDE 'TAU_COMMON' C C C Extensions to FORTRAN77 C ----------------------- C END DO / IMPLICIT NONE / INCLUDE / Names > 6 characters / C DO WHILE / VARIABLE FORMAT STATEMENTS / BYTE data type C C C Possible future upgrades C ------------------------ C C C Author/s C -------- C Jim Lewis RGO (jrl@uk.ac.cam.ast.mail) C C C History C ------- C 14-SEP-90 - Original program C 09-JAN-91 - Added in batch capabilities (JRL) C 12-SEP-91 - A few minor bug fixes. (JRL) C 06-JAN-92 - Fixed calls to TAUCAL routine. (JRL) C 28-JAN-92 - Minor modifications made for new 2d plotting routines. (JRL) C 01-FEB-94 - Modified for use with standalone UNIX FIGARO v4.0 (JRL) C C C Other Comments C -------------- C C+------------------------------------------------------------------------------ implicit none C C Local variables C integer address ! Dynamic memory address character access*6 ! Type of access for output file real alam0 ! Wavelenth zeropoint integer axslot ! Dyanmic memory axis slot logical badpix ! Flag for magic values logical badpixo ! Flag for magic values in output file logical batch_job ! TRUE if this is a batch job integer bytes_float ! Number of bytes in FLOAT type integer bytes_int ! Number of bytes in INT type integer dims(3) ! Dimensions of input cube real dlam ! Wavelength increment integer dptr ! Dynamic memory pointer for output data integer dslot ! Dyanmic memory slot for output data integer eptr ! Dynamic memory pointer for output error array integer eslot ! Dynamic memory slot for output error array real fsr ! Free spectral range integer i ! Integer counter integer iptr ! Dynamic memory pointer for input data integer islot ! Dynamic memory slot for input data integer ixptr ! Dynamic memory pointer for input X axis integer iyptr ! Dynamic memory pointer for input Y axis integer izptr ! Dynamic memory pointer for input Z axis byte magic_float ! Floating point magic value character name_fit*80 ! Name of output file character name_obs*80 ! Name of obs cube file integer ndim ! Number of dimension of input file integer nelm ! Number of elements in input file logical new_file ! Flag that output file is a new file. integer odims(3) ! Dimensions of output file character ref*5 ! Reference name of input cube integer status ! DSA status flag integer totspace ! Total space for workspace allocations character type*8 ! Data type of input file character typeo*8 ! Data type of output file integer wpt1 ! Dynamic memory pointer for workspace integer wpt2 ! Dynamic memory pointer for workspace integer wpt3 ! Dynamic memory pointer for workspace integer wslot ! Dynamic memory slot for workspace C C Maximum number of coefficients for a fit C integer ncomax parameter (ncomax=4) C C Functions C integer dsa_typesize integer dyn_element logical par_batch C C Common block with FSR C common /free/ fsr C C Numeric ranges. C include 'NDP_NUMERIC_RANGES' C C Common block for dynamic memory. C include 'DYNAMIC_MEMORY' C C TAUCAL common block C include 'TAU_COMMON' C C Open DSA system. C status = 0 call dsa_open(status) if (status .ne. 0) go to 500 C C Is this a batch job? C batch_job = par_batch() C C Open IMAGE file for calibration cube. C call par_rdchar('name_obs',' ',name_obs) call dsa_named_input('image',name_obs,status) call dsa_use_flagged_values('image',status) if (status .ne. 0) go to 500 C C Check image dimensions C call ndp_get_image_info('image',.true.,.false.,type,badpix,status) call dsa_data_size('image',3,ndim,dims,nelm,status) if (status .ne. 0) go to 500 if (ndim .ne. 3) then call dsa_wruser('Only 3-d data is acceptable...') call dsa_wrflush go to 500 end if C C Map IMAGE data array. C ref = 'image' call dsa_map_data('image','read','float',address,islot,status) if (status .ne. 0) go to 500 iptr = dyn_element(address) C C Map axis arrays C call dsa_map_axis_data('image',1,'read','float',address,axslot, : status) izptr = dyn_element(address) call dsa_map_axis_data('image',2,'read','float',address,axslot, : status) ixptr = dyn_element(address) call dsa_map_axis_data('image',3,'read','float',address,axslot, : status) iyptr = dyn_element(address) if (status .ne. 0) go to 500 C C Get parameters C call tau_cal_dta_name('image',name_obs,dta_name,lendta,status) call tau_cal_rdnum('image',0,'more.taurus.alam0','alam0',.true., : .true.,min_float,max_float,'angstroms',1,alam0,status) call tau_cal_rdnum('image',0,'more.taurus.dlam','dlam',.true., : .true.,min_float,max_float,'angstroms/frame',1,dlam,status) call tau_cal_rdnum('image',0,'more.taurus.fsr','fsr',.true., : .true.,min_float,max_float,'angstroms',1,fsr,status) if (status .ne. 0) go to 500 fsr = fsr/dlam C C Get some workspace for the options routine to use C bytes_float = dsa_typesize('float',status) bytes_int = dsa_typesize('int',status) totspace = 4*dims(3)*dims(2)*bytes_int + 2*ncomax*bytes_float call dsa_get_workspace(totspace,address,wslot,status) wpt1 = dyn_element(address) wpt2 = wpt1 + 4*dims(2)*dims(3)*bytes_int wpt3 = wpt2 + ncomax*bytes_float C C Now open output file.... C call par_rdkey('new_file',.true.,new_file) if (new_file) then call par_rdchar('name_fit',' ',name_fit) call dsa_named_output('output',name_fit,'image',1,1,status) odims(1) = ncomax + 4 odims(2) = dims(2) odims(3) = dims(3) call dsa_simple_output('output','data','float',3,odims,status) do i = 2,3 call dsa_reshape_axis('output',i,'image',i,1,odims(i), : status) end do call dsa_set_flagged_values('output',.true.,status) badpixo = .true. access = 'write' else call par_rdchar('name_fit',' ',name_fit) call dsa_named_input('output',name_fit,status) call ndp_get_image_info('output',.true.,.false.,typeo,badpixo, : status) call dsa_data_size('output',3,ndim,odims,nelm,status) if (status .ne. 0) go to 500 if (ndim .ne. 3) then call dsa_wruser('Only 3-d output file is acceptable...') call dsa_wrflush go to 500 end if if ((dims(2) .ne. odims(2)) .or. (dims(3) .ne. odims(3))) then call dsa_wruser('Output file must have the same spacial') call dsa_wruser('dimensions as the input file') call dsa_wrflush go to 500 end if if (odims(1) .ne. ncomax+4) then call dsa_wruser('Output file has wrong number of planes') call dsa_wrflush go to 500 end if access = 'update' end if if (badpixo) call dsa_use_flagged_values('output',status) call dsa_map_data('output',access,'float',address,dslot,status) dptr = dyn_element(address) call dsa_map_errors('output',access,'float',address,eslot,status) eptr = dyn_element(address) if (status .ne. 0) go to 500 C C Initialise output arrays if this is a new file... C if (new_file) then totspace = dims(2)*dims(3)*(ncomax+4) call dsa_get_flag_value('FLOAT',magic_float,status) if (status .ne. 0) go to 500 call gen_cfill(1,totspace,magic_float,dynamic_mem(dptr)) call gen_fill(totspace*bytes_float,0,dynamic_mem(eptr)) end if C C If this is an old file, then get the monitor array from the last C three frames of the output file otherwise, just fill it with zeros. C if (.not. new_file) then call rd_mon(dims(2),dims(3),ncomax,dynamic_mem(dptr), : dynamic_mem(wpt1)) else call gen_fill(4*dims(3)*dims(2)*bytes_int,0,dynamic_mem(wpt1)) end if C C Now do a little work...(rather a lot, really) C if (.not. batch_job) then call options(dynamic_mem(iptr),dims(2),dims(3),dims(1),ref, : dynamic_mem(wpt1),dynamic_mem(wpt2),dynamic_mem(wpt3), : ncomax,alam0,dlam,fsr,dynamic_mem(izptr),dynamic_mem(ixptr), : dynamic_mem(iyptr),dynamic_mem(dptr),dynamic_mem(eptr)) else call batch_opt(dynamic_mem(iptr),dims(2),dims(3),dims(1), : ncomax,dynamic_mem(wpt1),dynamic_mem(ixptr),dynamic_mem(iyptr), : fsr,alam0,dlam,dynamic_mem(wpt2),dynamic_mem(wpt3), : dynamic_mem(dptr),dynamic_mem(eptr)) end if C C Transfer the monitor array to the output file... C call set_mon(dims(2),dims(3),ncomax,dynamic_mem(wpt1), : dynamic_mem(dptr)) C C Tidy and exit C 500 call dsa_close(status) end subroutine options(input,nx,ny,nz,ref_name,mon,resmin, : resmax,ncomax,alam0,dlam,fsr,zaxis,xaxis,yaxis, : params,errpar) C+-------------------------------------------------------------------------------------- C OPTIONS --- Routine to present user with the options menu. C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) INPUT: (real array INPUT(NZ,NX,NY)) Input data cube. C (>) NX: (integer) X dimension of data C (>) NY: (integer) Y dimension of data C (>) NZ: (integer) Z dimension of data C (>) REF_NAME: (character) Reference name of input data cube. C (W) MON: (integer array MON(4,NX,NY)) Flags for marking spectra to be used C (W) RESMIN: (real array RESMIN(NCOMAX)) Array of result lower bounds C (W) RESMAX: (real array RESMAX(NCOMAX)) Array of result upper bounds C (>) NCOMAX: (integer) Maximum number of coefs for a fit C (>) ALAM0: (real) Wavelength zeropoint C (>) DLAM: (real) Wavelength increment C (>) FSR: (real) Free spectral range C (>) ZAXIS: (real array ZAXIS(NZ)) Z axis data C (>) XAXIS: (real array XAXIS(NX)) X axis data C (>) YAXIS: (real array YAXIS(NY)) Y axis data C (<) PARAMS: (real array PARAMS(NCOMAX+4,NX,NY)) Output fit coeffients C (<) ERRPAR: (real array ERRPAR(NCOMAX+4,NX,NY)) Errors in output coefs. C C C AUTHOR: Jim Lewis -- RGO C DATE: 14-SEP-1990 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nx,ny,nz,ncomax integer mon(4,nx,ny) real input(nz,nx,ny),resmin(ncomax), : resmax(ncomax),dlam,fsr,zaxis(nz),xaxis(nx),yaxis(ny), : params(ncomax+4,nx,ny),errpar(ncomax+4,nx,ny),alam0 character ref_name*(*) C C Functions C integer dyn_element integer dsa_typesize integer ich_fold integer ich_len logical par_qnum logical par_quest C C Local variables C integer idum,nco,nkeep,totspace,address,slot1,slot2,status,unit, : dptr,zptr,cptr,eptr,wpt2,wpt3,wpt4,bytes_float,bytes_int,deptr, : boxsize,mptr,slot3,nmx,nmy,wpt1,iwt,how_guess,nsptr,view_par, : ntriv_min,niter,ignore,odims(2),vptr,vslot,uptr,fit_type, : dat_slot,i,j,ipad1,ipad2,ntoss,wpt5,slot4,cvslot, : cpt1,cpt2,cpt3,cpt4,cpt5 real xmin,xmax,ymin,ymax,dummy,thresh(2),dchi_max, : minf,maxf,x1,x2,y1,y2,yrange(2) logical mask,picture,fixed_limits,changed,alloc, : use_curs,scale_y character opt*1,message*80,blank*1,helpfile*80 C C Common block for dynamic memory. C include 'DYNAMIC_MEMORY' C C Numeric ranges. C include 'NDP_NUMERIC_RANGES' C C Common block with convergence parameters C common /conv/ dchi_max,ntriv_min C C Common block with Y scaling info for line plotting routines C common /scaling/ scale_y,yrange C C A few initialisations C xmin = 1.0 xmax = float(nx) ymin = 1.0 ymax = float(ny) picture = .false. mask = .false. blank = ' ' changed = .true. alloc = .false. scale_y = .false. C C Initialise the fitting function C call var_getnum('type_fit',0,0,dummy,status) if ((status .ne. 0) .or. (dummy .eq. 0.0)) then status = 0 fit_type = 1 else fit_type = nint(dummy) end if call type_fit(.false.,fit_type,nco) C C Initialise weight scheme C call var_getnum('wt_scheme',0,0,dummy,status) if ((status .ne. 0) .or. (dummy .eq. 0.0)) then status = 0 iwt = 1 else iwt = nint(dummy) end if C C Initialise guessing algorithm C call var_getnum('how_guess',0,0,dummy,status) if ((status .ne. 0) .or. (dummy .eq. 0.0)) then status = 0 how_guess = 1 else how_guess = nint(dummy) if (how_guess .gt. 4) how_guess = 1 end if C C Initialise a few other things which may be parameters... C Start with smoothing box size. C call var_getnum('boxsize',0,0,dummy,status) if ((status .ne. 0) .or. (dummy .eq. 0.0)) then status = 0 boxsize = 1 else boxsize = nint(dummy) end if C C Number of iterations which chi-squared changes by a C trivial fraction in order for the solution to be considered C to be converged. C call var_getnum('ntriv',0,0,dummy,status) if ((status .ne. 0) .or. (dummy .eq. 0.0)) then status = 0 ntriv_min = 3 else ntriv_min = nint(dummy) end if C C Definition of a trivial change in the fraction of chi-squared. C call var_getnum('dchi',0,0,dummy,status) if ((status .ne. 0) .or. (dummy .eq. 0.0)) then status = 0 dchi_max = 0.05 else dchi_max = dummy end if C C Maximum number of iterations for any spectrum C call var_getnum('niter',0,0,dummy,status) if ((status .ne. 0) .or. (dummy .eq. 0.0)) then status = 0 niter = 30 else niter = nint(dummy) end if C C Set up some very useful constants C fixed_limits = .false. bytes_int = dsa_typesize('int',status) bytes_float = dsa_typesize('float',status) opt = ' ' call var_setchr('opt',0,0,blank,ignore) C C Tell the user what the current setup is like... C call info(fit_type,iwt,how_guess,boxsize,niter,ntriv_min,dchi_max) C C Options loop C do while (opt .ne. 'Q') C C Read option and change it to upper case C call dsa_wruser(' ') call dsa_wrflush call par_cnpar('opt') call par_rdchar('opt',' ',opt) call var_setchr('opt',0,0,blank,ignore) idum = ich_fold(opt) C C Help menu C if (opt .eq. 'H') then call dsa_open_text_file('TAU_FITS_OPT','.TXT','old', : .false.,unit,helpfile,status) if (status .ne. 0) then status = 0 call dsa_wruser('Sorry, I can''t find the HELP file') call dsa_wrflush go to 500 end if call dsa_wruser : ('****The options available in the menu are:') call dsa_wrflush call dsa_wruser(' ') call dsa_wrflush read(unit,105,iostat=status) message 105 format(a) do while (status .eq. 0) call dsa_wruser(message(:ich_len(message))) call dsa_wrflush read(unit,105,iostat=status) message end do call dsa_wruser(' ') call dsa_wrflush status = 0 close(unit) call dsa_free_lu(unit,status) C C Option to change the type of fit C else if (opt .eq. 'T') then call type_fit(.true.,fit_type,nco) C C Option to do the actual fitting C else if (opt .eq. 'F') then C C If you haven't fixed any limits, then there isn't much reason to C do this, is there? C if (.not. fixed_limits) then call dsa_wruser(' ') call dsa_wrflush call dsa_wruser('You haven''t fixed any limits yet. ') call dsa_wruser('When you have, let me know, but don''t') call dsa_wruser(' bug me until then...') call dsa_wrflush call dsa_wruser(' ') call dsa_wrflush C C Do this bit if you have changed the limits of your data or the smoothing C size since the last you binned the data... C else if (how_guess .eq. 5) then call dsa_wruser('You need to reinitialise the') call dsa_wruser(' guessing algorithm before ') call dsa_wruser(' proceding!') call dsa_wrflush go to 500 end if if (changed) then if (alloc) then call dsa_free_workspace(dat_slot,status) if (status .ne. 0) then status = 0 go to 500 end if alloc = .false. end if C C Throw out any spectra which may be on the very edge of the C cube as it will be impossible to smooth them properly C ntoss = 0 ipad1 = boxsize/2 if (mod(boxsize,2) .eq. 0) then ipad2 = ipad1 - 1 else ipad2 = ipad1 end if do j = 1,ny do i = 1,nx if ((i .le. ipad1) .or. (i .gt. nx-ipad2) .or. : (j .le. ipad1) .or. (j .gt. ny-ipad2)) then if (mon(1,i,j) .gt. 0) then mon(1,i,j) = 0 ntoss = ntoss + 1 end if end if end do end do if (ntoss .gt. 0) then write(message,107) ntoss 107 format(i6,' spectra thrown away on very edge') call dsa_wruser(message(:ich_len(message))) call dsa_wrflush end if nkeep = nkeep - ntoss C C Get workspace for data preparation routine C totspace = nz*2*nkeep*bytes_float + nz*bytes_int + : 3*nkeep*bytes_int call dsa_get_workspace(totspace,address,dat_slot, : status) if (status .ne. 0) then status = 0 go to 500 end if dptr = dyn_element(address) deptr = dptr + bytes_float*nz*nkeep nsptr = deptr + bytes_float*nz*nkeep uptr = nsptr + nz*bytes_int alloc = .true. C C Grab exactly the data you want and no more... C call prepare_data(input,nx,ny,nz,mon,nkeep,boxsize, : dynamic_mem(dptr),dynamic_mem(deptr), : dynamic_mem(nsptr),iwt,dynamic_mem(uptr)) changed = .false. end if C C Get some more workspace for the fitting routine C totspace = nco*(2*nkeep + 2*nco)*bytes_float + : nco*bytes_int + nz*bytes_float call dsa_get_workspace(totspace,address,slot2,status) if (status .ne. 0) then status = 0 go to 500 end if wpt2 = dyn_element(address) wpt3 = wpt2 + nco*nco*bytes_float wpt4 = wpt3 + nco*nco*bytes_float cptr = wpt4 + nco*bytes_int eptr = cptr + nco*nkeep*bytes_float zptr = eptr + nco*nkeep*bytes_float C C Do the fits C call fitt(nx,ny,nz,nkeep,mon,dynamic_mem(uptr), : dynamic_mem(dptr),nco,ncomax,fit_type,how_guess,fsr, : alam0,dlam,niter,picture,params,dynamic_mem(deptr), : dynamic_mem(zptr),dynamic_mem(wpt2),dynamic_mem(wpt3), : dynamic_mem(wpt4),dynamic_mem(cptr),dynamic_mem(eptr)) C C Write out the fits and free the workspace C call output(dynamic_mem(cptr),dynamic_mem(eptr),nx,ny, : nco,ncomax,nkeep,mon,dynamic_mem(uptr),fit_type,boxsize, : params,errpar) call dsa_free_workspace(slot2,status) end if C C Option to limit the regions for consideration in fitting C else if (opt .eq. 'L') then C C Find the limits of the regions C call limits(xmin,xmax,ymin,ymax,nx,ny,mask,thresh,xaxis, : yaxis,picture,mptr) fixed_limits = .true. changed = .true. C C If not using a mask, then create a bit of dummy workspace for the C next big subroutine call... C if (mask) then nmx = nx nmy = ny else call dsa_get_work_array(1,'float',address,slot3,status) if (status .ne. 0) then status = 0 go to 500 end if mptr = dyn_element(address) nmx = 1 nmy = 1 end if C C Prepare the 'monitor array' which tells the various routine which C spectra are to be used, and where in the data array they can be found. C call prepare_mon(nx,ny,nz,xmin,xmax,ymin,ymax, : nkeep,mask,thresh,nmx,nmy,dynamic_mem(mptr),mon) C C Message to the user C write(message,100) nkeep 100 format(i8,' spectra to be considered') call dsa_wruser(message(:ich_len(message))) call dsa_wrflush if (.not. mask) call dsa_free_workspace(slot3,status) C C Option to do limits in display mode (Picture) C else if (opt .eq. 'P') then call two_d_plot(nx,ny,xaxis,yaxis,picture,status) if (status .ne. 0) then status = 0 picture = .false. end if C C Option create a file from one of the results (Xtract) C else if (opt .eq. 'X') then call par_cnpar('xtract_field') call dsa_output('xtract_field','xtract_field',' ',1,1, : status) odims(1) = nx odims(2) = ny call dsa_simple_output('xtract_field','d,a1,a2','float',2, : odims,status) call dsa_reshape_axis('xtract_field',1,ref_name,2,1,nx, : status) call dsa_reshape_axis('xtract_field',2,ref_name,3,1,ny, : status) call dsa_set_flagged_values('xtract_field',.true.,status) call dsa_use_flagged_values('xtract_field',status) call dsa_map_data('xtract_field','write','float',address, : vslot,status) if (status .ne. 0) then call dsa_wruser('Trouble opening output file') call dsa_wrflush status = 0 go to 500 end if vptr = dyn_element(address) call dsa_wruser(' ') call dsa_wrflush call dsa_wruser('Choose the parameter to extract') call dsa_wrflush call dsa_wruser('(1) Peak') call dsa_wrflush call dsa_wruser('(2) Central wavelength') call dsa_wrflush call dsa_wruser('(3) Dispersion') call dsa_wrflush call dsa_wruser('(4) Background') call dsa_wrflush call dsa_wruser('(5) Error in Peak') call dsa_wrflush call dsa_wruser('(6) Error in Central Wavelength') call dsa_wrflush call dsa_wruser('(7) Error in Dispersion') call dsa_wrflush call dsa_wruser('(8) Error in Background') call dsa_wrflush call dsa_wruser('(9) Relative error in Peak') call dsa_wrflush call dsa_wruser : ('(10) Relative error in Central wavelength') call dsa_wrflush call dsa_wruser('(11) Relative error in Sigma') call dsa_wrflush call dsa_wruser('(12) Relative error in Background') call dsa_wrflush call dsa_wruser('(13) Velocity field') call dsa_wrflush call dsa_wruser('(14) Dispersion in velocity units') call dsa_wrflush call dsa_wruser('(Anything Else) None of the above') call dsa_wrflush call par_cnpar('xtract_par') call par_rdval('xtract_par',float(min_int),float(max_int), : 1.0,' ',dummy) view_par = nint(dummy) if ((view_par .lt. 1) .or. (view_par .gt. 14)) then call dsa_wruser('You obviously weren''t serious...') call dsa_wrflush else if (view_par .eq. 13) then call velocity_field(nx,ny,nz,ncomax,params, : dynamic_mem(vptr),minf,maxf) else call fillin_field(nx,ny,ncomax,view_par,params,errpar, : dynamic_mem(vptr),minf,maxf) end if call dsa_set_range('xtract_field',minf,maxf,status) call dsa_post_process_flagged_values('xtract_field',status) call dsa_unmap(vslot,status) call dsa_close_structure('xtract_field',status) C C Option to bin up in spatial direction C else if (opt .eq. 'B') then call par_cnpar('boxsize') call par_rdval('boxsize',float(min_int),float(max_int),1.0, : ' ',dummy) boxsize = nint(dummy) changed = .true. C C Option to use a mask image to restrict the action in the spatial direction C else if (opt .eq. 'M') then call par_cnpar('use_mask') call par_rdkey('use_mask',.true.,mask) C C Option to change weighting scheme... C else if (opt .eq. 'W') then call weights(iwt) C C Option to change guessing algorithm C else if (opt .eq. 'G') then call find_guess(how_guess,.false.,.false.) C C Option to change convergence parameters C else if (opt .eq. 'C') then call conv_pars(ntriv_min,dchi_max,niter) C C Reject bulk fits (Zap) C else if (opt .eq. 'Z') then call reject(ncomax,nco,nx,ny,picture,xaxis,yaxis,params, : errpar,mon,resmin,resmax) C C Edit individual fits C else if (opt .eq. 'E') then C C If you haven't fixed any limits, then there isn't much reason to C do this, is there? C if (.not. fixed_limits) then call dsa_wruser(' ') call dsa_wrflush call dsa_wruser('You haven''t fixed any limits yet. ') call dsa_wruser('When you have, let me know, but don''t') call dsa_wruser(' bug me until then...') call dsa_wrflush call dsa_wruser(' ') call dsa_wrflush else C C Do this bit if you have changed the limits of your data or the smoothing C size since the last you binned the data... C if (changed) then if (alloc) then call dsa_free_workspace(dat_slot,status) if (status .ne. 0) then status = 0 go to 500 end if alloc = .false. end if C C Throw out any spectra which may be on the very edge of the C cube as it will be impossible to smooth them properly C ntoss = 0 ipad1 = boxsize/2 if (mod(boxsize,2) .eq. 0) then ipad2 = ipad1 - 1 else ipad2 = ipad1 end if do j = 1,ny do i = 1,nx if ((i .le. ipad1) .or. (i .gt. nx-ipad2) .or. : (j .le. ipad1) .or. (j .gt. ny-ipad2)) then if (mon(1,i,j) .gt. 0) then mon(1,i,j) = 0 ntoss = ntoss + 1 end if end if end do end do if (ntoss .gt. 0) then write(message,107) ntoss call dsa_wruser(message(:ich_len(message))) end if nkeep = nkeep - ntoss C C Get some workspace for use in the data preparation routine C and the plotting routine C totspace = nz*2*nkeep*bytes_float + nz*bytes_int + : 3*nkeep*bytes_int call dsa_get_workspace(totspace,address,dat_slot, : status) if (status .ne. 0) then status = 0 go to 500 end if dptr = dyn_element(address) deptr = dptr + nz*nkeep*bytes_float nsptr = deptr + nz*nkeep*bytes_float uptr = nsptr + bytes_int*nz alloc = .true. C C Grab the data that you want. C call prepare_data(input,nx,ny,nz,mon,nkeep,boxsize, : dynamic_mem(dptr),dynamic_mem(deptr), : dynamic_mem(nsptr),iwt,dynamic_mem(uptr)) changed = .false. end if C C Get some more workspace for fitting routines C totspace = 3*nz*bytes_float call dsa_get_workspace(totspace,address,slot1,status) if (status .ne. 0) then status = 0 go to 500 end if zptr = dyn_element(address) wpt1 = zptr + bytes_float*nz wpt2 = wpt1 + bytes_float*nz C C Now edit the fits... C call edit(nx,ny,nz,nco,ncomax,picture, : dynamic_mem(zptr),params,errpar,dynamic_mem(dptr), : dynamic_mem(deptr),niter,nkeep,fit_type,how_guess, : alam0,dlam,mon,dynamic_mem(uptr),xmin,xmax,ymin,ymax, : dynamic_mem(wpt1),dynamic_mem(wpt2),iwt,boxsize) call dsa_free_workspace(slot1,status) end if C C Option for doing nearest neighbour fit... C else if (opt .eq. 'N') then C C If you haven't fixed any limits, then there isn't much reason to C do this, is there? C if (.not. fixed_limits) then call dsa_wruser(' ') call dsa_wrflush call dsa_wruser('You haven''t fixed any limits yet. ') call dsa_wruser('When you have, let me know, but don''t') call dsa_wruser(' bug me until then...') call dsa_wrflush call dsa_wruser(' ') call dsa_wrflush else C C Do this bit if you have changed the limits of your data or the smoothing C size since the last you binned the data... C if (changed) then if (alloc) then call dsa_free_workspace(dat_slot,status) if (status .ne. 0) then status = 0 go to 500 end if alloc = .false. end if C C Throw out any spectra which may be on the very edge of the C cube as it will be impossible to smooth them properly C ntoss = 0 ipad1 = boxsize/2 if (mod(boxsize,2) .eq. 0) then ipad2 = ipad1 - 1 else ipad2 = ipad1 end if do j = 1,ny do i = 1,nx if ((i .le. ipad1) .or. (i .gt. nx-ipad2) .or. : (j .le. ipad1) .or. (j .gt. ny-ipad2)) then if (mon(1,i,j) .gt. 0) then mon(1,i,j) = 0 ntoss = ntoss + 1 end if end if end do end do if (ntoss .gt. 0) then write(message,107) ntoss call dsa_wruser(message(:ich_len(message))) end if nkeep = nkeep - ntoss C C Get workspace for data preparation routine C totspace = nz*2*nkeep*bytes_float + nz*bytes_int + : 3*nkeep*bytes_int call dsa_get_workspace(totspace,address,dat_slot, : status) if (status .ne. 0) then status = 0 go to 500 end if dptr = dyn_element(address) deptr = dptr + bytes_float*nz*nkeep nsptr = deptr + bytes_float*nz*nkeep uptr = nsptr + nz*bytes_int alloc = .true. C C Grab exactly the data you want and no more... C call prepare_data(input,nx,ny,nz,mon,nkeep,boxsize, : dynamic_mem(dptr),dynamic_mem(deptr), : dynamic_mem(nsptr),iwt,dynamic_mem(uptr)) changed = .false. end if C C Grab a bit more workspace for fitting and plotting routines... C totspace = 2*nco*bytes_float + 4*nkeep*bytes_int call dsa_get_workspace(totspace,address,slot1,status) if (status .ne. 0) then status = 0 go to 500 end if cptr = dyn_element(address) eptr = cptr + nco*bytes_float wpt1 = eptr + nco*bytes_float wpt2 = wpt1 + 3*nkeep*bytes_int C C Now do the neighbour analysis C call neighbour(nx,ny,nz,ncomax,nkeep,fit_type,nco,niter, : boxsize,mon,dynamic_mem(uptr),params,dynamic_mem(dptr), : dynamic_mem(deptr),alam0,dlam,errpar,dynamic_mem(wpt1), : dynamic_mem(wpt2),resmin,resmax,dynamic_mem(cptr), : dynamic_mem(eptr),picture) call dsa_free_workspace(slot1,status) end if C C Option to display the results for a single spectrum (Show) C else if (opt .eq. 'S') then if (picture) then use_curs = par_quest('Use cursor to indicate limits?', : .true.) else use_curs = .false. end if if (use_curs) then call curse_limits(nx,ny,xaxis,yaxis,x1,x2,y1,y2,status) if (status .ne. 0) then status = 0 go to 500 end if else call dsa_wruser('Type coords for box') call dsa_wrflush ignore = par_qnum('X min',1.0,float(nx),1.0,.false., : ' ',x1) ignore = par_qnum('X max',1.0,float(nx),1.0,.false., : ' ',x2) ignore = par_qnum('Y min',1.0,float(ny),1.0,.false., : ' ',y1) ignore = par_qnum('Y max',1.0,float(ny),1.0,.false., : ' ',y2) if ((x1 .lt. 1.0) .or. (x2 .gt. float(nx)) .or. : (y1 .lt. 1.0) .or. (y2 .gt. float(ny))) go to 500 end if if (picture) then call pgwindow(1.0,float(nx),1.0,float(ny)) call pgmove(x1,y1) call pgdraw(x2,y1) call pgdraw(x2,y2) call pgdraw(x1,y2) call pgdraw(x1,y1) end if totspace = 4*ncomax*bytes_float call dsa_get_workspace(totspace,address,slot4,status) if (status .ne. 0) then status = 0 go to 500 end if wpt5 = dyn_element(address) call gen_fill(totspace,0,dynamic_mem(wpt5)) call show(nx,ny,x1,x2,y1,y2,params,errpar,ncomax,mon, : dynamic_mem(wpt5)) call dsa_free_workspace(slot4,status) if (status .ne. 0) then status = 0 go to 500 end if C C Option to display information about current setup (Info) C else if (opt .eq. 'I') then call info(fit_type,iwt,how_guess,boxsize,niter, : ntriv_min,dchi_max) C C Option to clean up velocity map (Velocities) C else if (opt .eq. 'V') then totspace = 2*nx*ny*bytes_float + 2*nx*ny*bytes_int + nx*ny call dsa_get_workspace(totspace,address,cvslot,status) if (status .ne. 0) then status = 0 go to 500 end if cpt1 = dyn_element(address) cpt2 = cpt1 + nx*ny cpt3 = cpt2 + nx*ny*bytes_float cpt4 = cpt3 + nx*ny*bytes_int cpt5 = cpt4 + nx*ny*bytes_float call cleanvel(nx,ny,ncomax,mon,params,xaxis,yaxis,fsr,alam0, : dlam,picture,ref_name,dynamic_mem(cpt1),dynamic_mem(cpt2), : dynamic_mem(cpt3),dynamic_mem(cpt4),dynamic_mem(cpt5)) call dsa_free_workspace(cvslot,status) if (status .ne. 0) then status = 0 go to 500 end if C C Option to put all line plots on the same scale (YSCALE) C else if (opt .eq. 'Y') then call par_cnpar('scale_y') call par_rdkey('scale_y',.false.,scale_y) if (scale_y) then call par_cnpar('yrange') call par_rdary('yrange',min_float,max_float,'a',' ', : 2,2,yrange) end if end if 500 end do end subroutine batch_opt(input,nx,ny,nz,ncomax,mon,xaxis,yaxis,fsr, : alam0,dlam,resmin,resmax,params,errpar) C+-------------------------------------------------------------------------------------- C BATCH_OPT --- Routine to do the necessary preparations while in batch C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) INPUT: (real array INPUT(NZ,NX,NY)) Input data cube. C (>) NX: (integer) X dimension of data C (>) NY: (integer) Y dimension of data C (>) NZ: (integer) Z dimension of data C (>) NCOMAX: (integer) Maximum number of coefs for a fit C (W) MON: (integer array MON(4,NX,NY)) Flags for marking spectra to be used C (>) XAXIS: (real array XAXIS(NX)) X axis data C (>) YAXIS: (real array YAXIS(NY)) Y axis data C (>) FSR: (real) Free spectral range C (>) ALAM0: (real) Wavelength zeropoint C (>) DLAM: (real) Wavelength increment C (W) RESMIN: (real array RESMIN(NCOMAX)) Array of result lower bounds C (W) RESMAX: (real array RESMAX(NCOMAX)) Array of result upper bounds C (<) PARAMS: (real array PARAMS(NCOMAX+4,NX,NY)) Output fit coeffients C (<) ERRPAR: (real array ERRPAR(NCOMAX+4,NX,NY)) Errors in output coefs. C C C AUTHOR: Jim Lewis -- RGO C DATE: 09-JAN-1991 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nx,ny,nz,ncomax integer mon(4,nx,ny) real xaxis(nx),yaxis(ny),fsr,alam0,dlam,params(ncomax+4,nx,ny), : errpar(ncomax+4,nx,ny),resmin(ncomax),resmax(ncomax), : input(nz,nx,ny) C C Functions C integer dsa_typesize integer dyn_element integer ich_len logical par_given C C Local variables C integer boxsize,ntriv_min,niter,fit_type,nco,iwt,mptr,nmx,nmy, : address,slot,status,nkeep,batch_fit,how_guess,ntoss,ipad1,ipad2, : i,j,bytes_float,bytes_int,totspace,dat_slot,dptr,deptr,nsptr, : uptr,wpt1,wpt2,wpt3,wpt4,cptr,eptr,zptr,nsweep real dummy,dchi_max,xmin,xmax,ymin,ymax,thresh(2),dum(2) logical given,mask,given2,picture character message*80 C C Parameters for DSA and PAR routines... C character range(4)*7 data range /'rng_amp','rng_pos','rng_sig','rng_bac'/ C C Common block for dynamic memory. C include 'DYNAMIC_MEMORY' C C Numeric ranges. C include 'NDP_NUMERIC_RANGES' C C Common block with convergence parameters C common /conv/ dchi_max,ntriv_min C C Begin by initialising some variables. If the variable hasn't been C given in the command line, then the program will try to find it in C the BVARS.DST file. If it's not there, then an error occurs and the C routine quits. C C First the smoothing box size C given = par_given('boxsize') if (given) then call par_rdval('boxsize',1.0,float(max_int),1.0,'pixels',dummy) boxsize = nint(dummy) else call var_getnum('boxsize',0,0,dummy,status) boxsize = nint(dummy) if ((status .ne. 0) .or. (boxsize .lt. 1)) then call dsa_wruser('BOXSIZE parameter undefined') call dsa_wrflush go to 500 end if end if C C Number of iterations which chi-squared changes by a C trivial fraction in order for the solution to be considered C to be converged. C given = par_given('ntriv') if (given) then call par_rdval('ntriv',1.0,float(max_int),2.0,'iterations', : dummy) ntriv_min = nint(dummy) else call var_getnum('ntriv',0,0,dummy,status) ntriv_min = nint(dummy) if ((status .ne. 0) .or. (ntriv_min .lt. 1)) then call dsa_wruser('NTRIV parameter undefined') call dsa_wrflush go to 500 end if end if C C Definition of a trivial change in the fraction of chi-squared. C given = par_given('dchi') if (given) then call par_rdval('dchi',0.0,max_float,0.05,' ',dchi_max) else call var_getnum('dchi',0,0,dchi_max,status) if ((status .ne. 0) .or. (dchi_max .le. 0.0)) then call dsa_wruser('DCHI parameter undefined') call dsa_wrflush go to 500 end if end if C C Maximum number of iterations for any spectrum C given = par_given('niter') if (given) then call par_rdval('niter',1.0,float(max_int),30.0,'iterations', : dummy) niter = nint(dummy) else call var_getnum('niter',0,0,dummy,status) niter = nint(dummy) if ((status .ne. 0) .or. (niter .lt. 1)) then call dsa_wruser('NITER parameter undefined') call dsa_wrflush go to 500 end if end if C C Initialise the fitting function C given = par_given('type_fit') if (given) then call par_rdval('type_fit',1.0,2.0,1.0,' ',dummy) fit_type = nint(dummy) else call var_getnum('type_fit',0,0,dummy,status) fit_type = nint(dummy) if ((status .ne. 0) .or. (fit_type .lt. 1) .or. : (fit_type .gt. 2)) then call dsa_wruser('TYPE_FIT parameter undefined') call dsa_wrflush go to 500 end if end if call type_fit(.false.,fit_type,nco) C C Initialise weight scheme C given = par_given('wt_scheme') if (given) then call par_rdval('wt_scheme',1.0,2.0,1.0,' ',dummy) iwt = nint(dummy) else call var_getnum('wt_scheme',0,0,dummy,status) iwt = nint(dummy) if ((status .ne. 0) .or. (iwt .lt. 1) .or. (iwt .gt. 2)) then call dsa_wruser('WT_SCHEME parameter undefined') call dsa_wrflush go to 500 end if end if C C Write out all this info now C call info(fit_type,iwt,how_guess,boxsize,niter,ntriv_min,dchi_max) C C Check to see if a mask is being used. If it is not specified on C the command line, then assume it isn't being used C given = par_given('use_mask') if (given) then call par_rdkey('use_mask',.true.,mask) if (mask) then given2 = par_given('name_mask') .and. : par_given('mask_thresh') if (.not. given2) then call dsa_wruser('Either NAME_MASK or MASK_THRESH') call dsa_wruser(' parameters undefined') call dsa_wrflush go to 500 end if end if else mask = .false. end if C C No pictures allowed in batch! C picture = .false. C C See if the x,y limits have been given on the command line. If C not, then quit, since you are unlikely to be able to allocate enough C workspace to do fits on the whole frame. C given2 = par_given('start') .and. par_given('end') if (.not. given2) then call dsa_wruser('START or END parameters undefined') call dsa_wrflush go to 500 end if C C Find the limits of the regions C call limits(xmin,xmax,ymin,ymax,nx,ny,mask,thresh,xaxis, : yaxis,picture,mptr) C C If not using a mask, then create a bit of dummy workspace for the C next big subroutine call... C if (mask) then nmx = nx nmy = ny else call dsa_get_work_array(1,'float',address,slot,status) if (status .ne. 0) then status = 0 go to 500 end if mptr = dyn_element(address) nmx = 1 nmy = 1 end if C C Prepare the 'monitor array' which tells the various routine which C spectra are to be used, and where in the data array they can be found. C call prepare_mon(nx,ny,nz,xmin,xmax,ymin,ymax, : nkeep,mask,thresh,nmx,nmy,dynamic_mem(mptr),mon) C C Message to the user C write(message,100) nkeep 100 format(i8,' spectra to be considered') call dsa_wruser(message(:ich_len(message))) call dsa_wrflush if (.not. mask) call dsa_free_workspace(slot,status) C C Find out which you want to do, (1) fit, with some sort of guessing or C (2) nearest neighbour fitting. C given = par_given('batch_fit') if (.not. given) then call dsa_wruser : ('BATCH_FIT parameter unspecified on command line') call dsa_wrflush go to 500 else call par_rdval('batch_fit',1.0,2.0,1.0,' ',dummy) batch_fit = nint(dummy) end if C C If doing nearest neighbour, then check to see if various parameters C are either given on the command line or are included in BVARS.DST C if (batch_fit .eq. 2) then given = par_given('nsweep') if (.not. given) then call var_getnum('nsweep',0,0,dummy,status) nsweep = nint(dummy) if ((status .ne. 0) .or. (nsweep .le. 0)) then call dsa_wruser('NSWEEP parameter undefined') call dsa_wrflush go to 500 end if end if do i = 1,nco given = par_given(range(i)) if (.not. given) then call var_getnum(range(i),1,2,dum,status) if (status .ne. 0) then call dsa_wruser(range(i)//' parameters undefined') call dsa_wrflush go to 500 end if end if end do C C Otherwise, check to make sure you have included a guessing algorithm C else given = par_given('how_guess') if (.not. given) then call var_getnum('how_guess',1,2,dum,status) if (status .ne. 0) then call dsa_wruser('HOW_GUESS parameter undefined') call dsa_wrflush go to 500 end if end if end if C C Prepare the data which you want to fit. First throw out any C spectra which may be on the very edge of the cube as it will C be impossible to smooth them properly C ntoss = 0 ipad1 = boxsize/2 if (mod(boxsize,2) .eq. 0) then ipad2 = ipad1 - 1 else ipad2 = ipad1 end if do j = 1,ny do i = 1,nx if ((i .le. ipad1) .or. (i .gt. nx-ipad2) .or. : (j .le. ipad1) .or. (j .gt. ny-ipad2)) then if (mon(1,i,j) .gt. 0) then mon(1,i,j) = 0 ntoss = ntoss + 1 end if end if end do end do if (ntoss .gt. 0) then write(message,107) ntoss 107 format(i6,' spectra thrown away on very edge') call dsa_wruser(message(:ich_len(message))) call dsa_wrflush end if nkeep = nkeep - ntoss C C Get workspace for data preparation routine C bytes_float = dsa_typesize('float',status) bytes_int = dsa_typesize('int',status) totspace = nz*2*nkeep*bytes_float + nz*bytes_int + : 3*nkeep*bytes_int call dsa_get_workspace(totspace,address,dat_slot,status) if (status .ne. 0) then call dsa_wruser : ('Couldn''t allocate workspace for data preparation') call dsa_wrflush go to 500 end if dptr = dyn_element(address) deptr = dptr + bytes_float*nz*nkeep nsptr = deptr + bytes_float*nz*nkeep uptr = nsptr + nz*bytes_int C C Grab exactly the data you want and no more... C call prepare_data(input,nx,ny,nz,mon,nkeep,boxsize, : dynamic_mem(dptr),dynamic_mem(deptr),dynamic_mem(nsptr),iwt, : dynamic_mem(uptr)) C C If you are doing a fit with guessing then get the guessing algorithm C if (batch_fit .eq. 1) then call find_guess(how_guess,.true.,.false.) C C Get some more workspace for the fitting routine C totspace = nco*(2*nkeep + 2*nco)*bytes_float + : nco*bytes_int + nz*bytes_float call dsa_get_workspace(totspace,address,slot,status) if (status .ne. 0) then call dsa_wruser : ('Couldn''t allocate workspace for fitting') call dsa_wrflush go to 500 end if wpt2 = dyn_element(address) wpt3 = wpt2 + nco*nco*bytes_float wpt4 = wpt3 + nco*nco*bytes_float cptr = wpt4 + nco*bytes_int eptr = cptr + nco*nkeep*bytes_float zptr = eptr + nco*nkeep*bytes_float C C Do the fits C call fitt(nx,ny,nz,nkeep,mon,dynamic_mem(uptr), : dynamic_mem(dptr),nco,ncomax,fit_type,how_guess,fsr, : alam0,dlam,niter,picture,params,dynamic_mem(deptr), : dynamic_mem(zptr),dynamic_mem(wpt2),dynamic_mem(wpt3), : dynamic_mem(wpt4),dynamic_mem(cptr),dynamic_mem(eptr)) C C Write out the fits and free the workspace C call output(dynamic_mem(cptr),dynamic_mem(eptr),nx,ny, : nco,ncomax,nkeep,mon,dynamic_mem(uptr),fit_type,boxsize,params, : errpar) C C If you are doing a nearest neighbour fit, then start off by getting C a bit of workspace for the fitting routine... C else if (batch_fit .eq. 2) then totspace = 2*nco*bytes_float + 4*nkeep*bytes_int call dsa_get_workspace(totspace,address,slot,status) if (status .ne. 0) then status = 0 go to 500 end if cptr = dyn_element(address) eptr = cptr + nco*bytes_float wpt1 = eptr + nco*bytes_float wpt2 = wpt1 + 3*nkeep*bytes_int C C Now do the neighbour analysis C call neighbour(nx,ny,nz,ncomax,nkeep,fit_type,nco,niter, : boxsize,mon,dynamic_mem(uptr),params,dynamic_mem(dptr), : dynamic_mem(deptr),alam0,dlam,errpar,dynamic_mem(wpt1), : dynamic_mem(wpt2),resmin,resmax,dynamic_mem(cptr), : dynamic_mem(eptr),picture) end if C C Free the workspaces and exit C call dsa_free_workspace(dat_slot,status) call dsa_free_workspace(slot,status) 500 end subroutine show(nx,ny,x1,x2,y1,y2,params,errpar,ncomax,mon, : work) C+-------------------------------------------------------------------------------------- C SHOW --- Show a set of fit parameters to the user C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) NX: (integer) X dimension of frames C (>) NY: (integer) Y dimension of frames C (>) X1: (integer) Min X coord of region in question C (>) X2: (integer) Max X coord of region in question C (>) Y1: (integer) Min Y coord of region in question C (>) Y2: (integer) Max Y coord of region in question C (>) PARAMS: (real array PARAMS(NCOMAX+4,NX,NY)) Output data array C (>) ERRPAR: (real array ERRPAR(NCOMAX+4,NX,NY)) Error in output array C (>) NCOMAX: (integer) Maximum number of coefs for fits C (>) MON: (integer array MON(4,NX,NY)) Monitor array C (W) WORK: (real array WORK(4,NCOMAX)) Workspace C C C AUTHOR: Jim Lewis -- RGO C DATE: 19-DEC-1990 C UPDATE: 24-FEB-1991 -- Does stats over a box rather than for single C spectra. (JRL) C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nx,ny,ncomax integer mon(4,nx,ny) real params(ncomax+4,nx,ny),errpar(ncomax+4,nx,ny),work(4,ncomax), : x1,x2,y1,y2 C C Functions C integer ich_len C C Local variables C integer fit_type,nco,idone,i,ngood,ix,iy real val character message*80 C C Numeric ranges C include 'NDP_NUMERIC_RANGES' C C Initialise part of the work array C do i = 1,ncomax work(3,i) = max_float work(4,i) = min_float end do C C Sum up all of the good results in this region C ngood = 0 do iy = nint(y1),nint(y2) do ix = nint(x1),nint(x2) C C First find out if there has been a successful fit... C idone = mon(3,ix,iy) C C If the fit has been done successfully, then add it it C if (idone .gt. 0) then fit_type = mon(2,ix,iy) call type_fit(.false.,fit_type,nco) ngood = ngood + 1 do i = 1,nco val = params(i,ix,iy) work(1,i) = work(1,i) + val work(2,i) = work(2,i) + errpar(i,ix,iy) work(3,i) = min(work(3,i),val) work(4,i) = max(work(4,i),val) end do end if end do end do C C Write out results. Check first to make sure there is at least one good one! C if (ngood .lt. 1) then call dsa_wruser('None of the spectra in this region have') call dsa_wruser(' been successfully fit!') call dsa_wrflush else do i = 1,ncomax work(1,i) = work(1,i)/float(ngood) work(2,i) = work(2,i)/float(ngood) end do write(message,100) (work(3,i),i=1,ncomax) 100 format('Mins: ',f12.5) call dsa_wruser(message(:ich_len(message))) call dsa_wrflush write(message,101) (work(4,i),i=1,ncomax) 101 format('Maxs: ',f12.5) call dsa_wruser(message(:ich_len(message))) call dsa_wrflush write(message,102) (work(1,i),i=1,ncomax) 102 format(':',f12.5) call dsa_wruser(message(:ich_len(message))) call dsa_wrflush write(message,103) (work(2,i),i=1,ncomax) 103 format(': ',f12.5) call dsa_wruser(message(:ich_len(message))) call dsa_wrflush end if end subroutine rd_mon(nx,ny,ncomax,params,mon) C+-------------------------------------------------------------------------------------- C RD_MON --- Routine to read the monitor array from an old C output file C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) NX: (integer) X dimension of frames C (>) NY: (integer) Y dimension of frames C (>) NCOMAX: (integer) Maximum number of coefs for fits C (>) PARAMS: (real array PARAMS(NCOMAX+4,NX,NY)) Output data array C (<) MON: (integer array MON(4,NX,NY)) Output monitor array C C C AUTHOR: Jim Lewis -- RGO C DATE: 27-NOV-1990 C UPDATE: 01-FEB-1994 -- Replaced call to cnv_fmtcon with call to C dsa_fmtcon (JRL) C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nx,ny,ncomax integer mon(4,nx,ny) real params(ncomax+4,nx,ny) C C Local variables C integer i,j,k,val_i,nbad real val_r C C Transfer the last three frames of output and convert to integer C do j = 1,ny do i = 1,nx do k = 1,4 val_r = params(k+ncomax,i,j) call dsa_fmtcon(.true.,4,3,val_r,val_i,1,nbad) mon(k,i,j) = val_i end do end do end do end subroutine set_mon(nx,ny,ncomax,mon,params) C+-------------------------------------------------------------------------------------- C SET_MON --- Routine to write the monitor array to the C output file C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) NX: (integer) X dimension of frames C (>) NY: (integer) Y dimension of frames C (>) NCOMAX: (integer) Maximum number of coefs for fits C (>) MON: (integer array MON(4,NX,NY)) Output monitor array C (!) PARAMS: (real array PARAMS(NCOMAX+4,NX,NY)) Output data array C C C AUTHOR: Jim Lewis -- RGO C DATE: 27-NOV-1990 C UPDATE: 01-FEB-1994 -- Replaced call to cnv_fmtcon with call to C dsa_fmtcon (JRL) C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nx,ny,ncomax integer mon(4,nx,ny) real params(ncomax+4,nx,ny) C C Local variables C integer i,j,nbad C C Transfer monitor array to output and convert to real C do j = 1,ny do i = 1,nx call dsa_fmtcon(.true.,3,4,mon(1,i,j), : params(ncomax+1,i,j),4,nbad) end do end do end subroutine fillin_field(nx,ny,ncomax,nfield,params,errpar, : field,amin,amax) C+-------------------------------------------------------------------------------------- C FILLIN_FIELD --- Routine to write results to a 2-d array C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) NX: (integer) X dimension of frames C (>) NY: (integer) Y dimension of frames C (>) NCOMAX: (integer) Maximum number of coefs for fits C (>) NFIELD: (integer) Number of the parameter desired C (>) PARAMS: (real array PARAMS(NCOMAX+4,NX,NY)) Output data array C (>) ERRPAR (real array ERRPAR(NCOMAX+4,NX,NY)) Output error array C (<) FIELD: (real array FIELD(NX,NY)) Output data array C (<) AMIN: (real) Minimum of output field (excluding magic values) C (<) AMAX: (real) Maximum of output field (excluding magic values) C C C AUTHOR: Jim Lewis -- RGO C DATE: 27-NOV-1990 C UPDATE: 19-DEC-1990 -- Now passes min and max back. (JRL) C 24-FEB-1991 -- Does sigma as a velocity now. (JRL) C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nx,ny,ncomax,nfield real params(ncomax+4,nx,ny),field(nx,ny),errpar(ncomax+4,nx,ny), : amin,amax C C Functions C integer ich_len C C Local variables C integer i,j,status real val,val1,val2,c,magic_float character message*80 C C Numeric ranges C include 'NDP_NUMERIC_RANGES' C C Initialise min and max C amin = max_float amax = min_float c = 3.0e5 call dsa_get_flag_value('FLOAT',magic_float,status) if (status .ne. 0) go to 500 C C Transfer the data and find the min and max C do j = 1,ny do i = 1,nx if (nfield .le. 4) then val = params(nfield,i,j) else if (nfield .le. 8) then val = errpar(nfield-4,i,j) else if (nfield .le. 12) then val1 = errpar(nfield-8,i,j) val2 = params(nfield-8,i,j) if ((val1 .gt. magic_float) .and. : (val2 .gt. magic_float)) then val = abs(val1/val2) else val = magic_float end if else if (nfield .eq. 14) then val1 = params(3,i,j) val2 = params(2,i,j) if ((val1 .gt. magic_float) .and. : (val2 .gt. magic_float)) then val = (val1/val2)*c else val = magic_float end if end if field(i,j) = val if (val .gt. magic_float) then amin = min(amin,val) amax = max(amax,val) end if end do end do C C Write out the min and max C write(message,100) amin,amax 100 format('The range is:',2(1pg12.4)) call dsa_wrflush call dsa_wruser(message(:ich_len(message))) 500 end subroutine velocity_field(nx,ny,nz,ncomax,params,result,min_vel, : max_vel) C+-------------------------------------------------------------------------------------- C VELOCITY_FIELD --- Routine to write velocity results to a 2-d array C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) NX: (integer) X dimension of frames C (>) NY: (integer) Y dimension of frames C (>) NZ: (integer) Number of pixels in spectral direction C (>) NCOMAX: (integer) Maximum number of coefs for fits C (>) PARAMS: (real array PARAMS(NCOMAX+4,NX,NY)) Output data array C (<) RESULT: (real array RESULT(NX,NY)) Output data array C (<) MIN_VEL: (real) Minimum velocity (excluding magic values) C (<) MAX_VEL: (real) Maximum velocity (excluding magic values) C C C AUTHOR: Jim Lewis -- RGO C DATE: 27-NOV-1990 C UPDATE: 19-DEC-1990 -- Now passes back min and max. Also asks for C rest wavelength. (JRL) C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nx,ny,nz,ncomax real params(ncomax+4,nx,ny),result(nx,ny),min_vel,max_vel C C Functions C integer ich_len logical par_batch C C Local variables C integer i,j,status real lam_cen,c,lamda,vel,magic_float character message*80 C C Numeric ranges C include 'NDP_NUMERIC_RANGES' C C Get central wavelength. Set speed of light (km/sec) C if (.not. par_batch()) call par_cnpar('lam_cen') call par_rdval('lam_cen',min_float,max_float,1.0,'angstroms', : lam_cen) c = 3.0e5 min_vel = max_float max_vel = min_float call dsa_get_flag_value('FLOAT',magic_float,status) if (status .ne. 0) go to 500 C C Now do the calculations C do j = 1,ny do i = 1,nx lamda = params(2,i,j) if (lamda .gt. magic_float) then vel = c*(lamda - lam_cen)/lam_cen result(i,j) = vel min_vel = min(vel,min_vel) max_vel = max(vel,max_vel) else result(i,j) = magic_float end if end do end do C C Give the user some indication of the range of velocities C with which he/she/it has to deal... C write(message,100) min_vel,max_vel 100 format('The range in velocity is:',2(1pg12.4),' km/sec') call dsa_wruser(message(:ich_len(message))) call dsa_wrflush 500 end subroutine conv_pars(ntriv_min,dchi_max,niter) C+-------------------------------------------------------------------------------------- C CONV_PARS --- Routine to set parameters for convergence of fits C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (!) NTRIV_MIN: (integer) Minimum number of trivial fractional changes to C chi-squared for convergence. C (!) DCHI_MAX: (integer) Definition of trivial fractional change C (!) NITER: (integer) Maximum number of iterations C C C AUTHOR: Jim Lewis -- RGO C DATE: 27-NOV-1990 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer ntriv_min,niter real dchi_max C C Functions C logical par_batch C C Local variables C real dummy logical batch C C Numeric ranges. C include 'NDP_NUMERIC_RANGES' C C Prompt for the convergence parameters C batch = par_batch() if (.not. batch) call par_cnpar('ntriv') call par_rdval('ntriv',1.0,float(max_int),1.0,' ',dummy) ntriv_min = nint(dummy) if (.not. batch) call par_cnpar('dchi') call par_rdval('dchi',0.0,max_float,1.0,' ',dchi_max) if (.not. batch) call par_cnpar('niter') call par_rdval('niter',0.0,max_float,1.0,' ',dummy) niter = nint(dummy) end subroutine find_guess(how_guess,batch,edit) C+-------------------------------------------------------------------------------------- C FIND_GUESS --- Routine to allow user to change the guessing algorithm C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C C (!) HOW_GUESS: (integer) Number of the guessing algorithm to be used C (>) BATCH: (logical) TRUE if this is a batch job C (>) EDIT: (logical) TRUE if this is called from EDIT. C C C AUTHOR: Jim Lewis -- RGO C DATE: 14-SEP-1990 C UPDATE: 09-JAN-1991 -- BATCH parameter added. (JRL) C 10-JAN-1991 -- EDIT parameter added (JRL) C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer how_guess logical batch,edit C C Local variables C real dummy,max C C Get the guessing scheme flag C call dsa_wruser(' ') call dsa_wrflush call dsa_wruser : ('You have a choice of the following guessing algorithms:') call dsa_wrflush call dsa_wruser('(1) Moments') call dsa_wrflush call dsa_wruser('(2) Previous Results') call dsa_wrflush if (batch) then max = 2.0 else call dsa_wruser('(3) By Hand (slow version)') call dsa_wrflush call dsa_wruser('(4) By Hand (fast version)') call dsa_wrflush max = 4.0 if (edit) then call dsa_wruser('(5) From another spectrum') call dsa_wrflush max = 5.0 end if end if if (.not. batch) call par_cnpar('how_guess') call par_rdval('how_guess',1.0,max,1.0,' ',dummy) how_guess = nint(dummy) end subroutine weights(iwt) C+-------------------------------------------------------------------------------------- C WEIGHTS --- Routine to allow user to change the way weights are assigned C to the data points C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C C (!) IWT: (integer) Number of the weighting scheme C C C AUTHOR: Jim Lewis -- RGO C DATE: 14-SEP-1990 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer iwt C C Functions C logical par_batch C C Local variables C real dummy C C Get weighting scheme flag C call dsa_wruser(' ') call dsa_wrflush call dsa_wruser : ('You may weight the data in the following ways:') call dsa_wrflush call dsa_wruser('(1) Sqrt(n)') call dsa_wrflush call dsa_wruser('(2) Equal weights') call dsa_wrflush if (.not. par_batch()) call par_cnpar('wt_scheme') call par_rdval('wt_scheme',1.0,2.0,1.0,' ',dummy) iwt = nint(dummy) end subroutine type_fit(ask,fit_type,nco) C+-------------------------------------------------------------------------------------- C TYPE_FIT --- Routine to allow user to change the fitting function used C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) ASK: (integer) TRUE if the user is to be prompted C (!) FIT_TYPE: (integer) Number of the fitting function C (<) NCO: (integer) Number of fit coefs for the fitting functions C C C AUTHOR: Jim Lewis -- RGO C DATE: 14-SEP-1990 C UPDATE: 19-DEC-1990 -- Added ASK parameter (JRL) C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nco,fit_type logical ask C C Functions C logical par_batch C C Local variables C integer ncoef(2) real dummy data ncoef /4,4/ C C Numeric ranges C include 'NDP_NUMERIC_RANGES' C C Select the type of fit and determine the number of fit C coefs for that fitting function C if (ask) then call dsa_wruser(' ') call dsa_wrflush call dsa_wruser('Which function to use in fits?') call dsa_wrflush call dsa_wruser('(1) Gaussian') call dsa_wrflush call dsa_wruser('(2) Cauchy') call dsa_wrflush if (.not. par_batch()) call par_cnpar('type_fit') call par_rdval('type_fit',1.0,2.0,1.0,' ',dummy) fit_type = nint(dummy) end if nco = ncoef(fit_type) end subroutine prepare_mon(nx,ny,nz,xmin,xmax,ymin,ymax,nkeep, : mask,thresh,nmx,nmy,mask_val,mon) C+-------------------------------------------------------------------------------------- C PREPARE_MON --- Routine to prepare the array which flags which spectra C will be used. This array will later also tell where in the C data workspace the spectrum can be found (but that is not C set up here...) C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) NX: (integer) X dimension of the data C (>) NY: (integer) Y dimension of the data C (>) NZ: (integer) Z dimension of the data C (>) XMIN: (real) Minimum value of X for consideration C (>) XMAX: (real) Maximum value of X for consideration C (>) YMIN: (real) Minimum value of Y for consideration C (>) YMAX: (real) Maximum value of Y for consideration C (<) NKEEP: (integer) Total number of spectra to be used C (>) MASK: (logical) Instruction to use a mask image C (>) THRESH: (real array THRESH(2)) Min and max mask thresholds C (>) NMX: (integer) X dimension of mask image C (>) NMY: (integer) Y dimension of mask image C (>) MASK_VAL: (real array MASK(NX,NY)) Input mask image C (<) MON: (integer array MON(4,NX,NY)) Output 'monitor' array C C C AUTHOR: Jim Lewis -- RGO C DATE: 14-SEP-1990 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nx,ny,nz,nkeep,nmx,nmy integer mon(4,nx,ny) real xmin,xmax,ymin,ymax,thresh(2),mask_val(nmx,nmy) logical mask C C Local variables C integer i,j real x,y logical keep C C First initialise monitor array C do i = 1,ny do j = 1,nx mon(1,j,i) = -1 end do end do C C Now create the monitor array C nkeep = 0 do i = 1,ny y = float(i) do j = 1,nx keep = .true. x = float(j) C C If you're using a mask see if the point lies between C the two mask thresholds. C if (mask) then if ((mask_val(j,i) .lt. thresh(1)) .or. : (mask_val(j,i) .gt. thresh(2))) keep = .false. end if C C See if the point is between the x,y limits C if ((x .gt. xmax) .or. (x .lt. xmin) .or. (y .gt. ymax) : .or. (y .lt. ymin)) keep = .false. C C Count the number of spectra kept and flag the ones that aren't C if (keep) then nkeep = nkeep + 1 else mon(1,j,i) = 0 end if end do end do end subroutine prepare_data(input,nx,ny,nz,mon,nkeep,boxsize,data, : err,nsum,iwt,used) C+-------------------------------------------------------------------------------------- C PREPARE_DATA --- Routine to prepare an array which contains only the data C which the user wants to use at a given time. C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) INPUT: (real array INPUT(NZ,NX,NY)) Input data cube. C (>) NX: (integer) X dimension of the data C (>) NY: (integer) Y dimension of the data C (>) NZ: (integer) Z dimension of the data C (>) MON: (integer array MON(4,NX,NY)) Input 'monitor' array C (>) NKEEP: (integer) Total number of spectra to be used C (>) BOXSIZE: (integer) Size of smoothing box. C (<) DATA: (real array DATA(NZ,NKEEP)) Output data array C (<) ERR: (real array ERR(NZ,NKEEP)) Output error array C (W) NSUM: (integer array NSUM(NZ)) Workspace C (>) IWT: (integer) Weighting scheme flag C (<) USED: (integer array USED(3,NKEEP)) Array with cell coords of used spectra C C C AUTHOR: Jim Lewis -- RGO C DATE: 14-SEP-1990 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nx,ny,nz,nkeep,boxsize,iwt integer mon(4,nx,ny),nsum(nz),used(3,nkeep) real input(nz,nx,ny),data(nz,nkeep),err(nz,nkeep) C C Functions C integer dsa_typesize C C Local variables C integer i,j,k,n,ii,ipad1,ipad2,totspace,jj,status,iz,nbad,status real d,idata,magic_float C C Numeric ranges C include 'NDP_NUMERIC_RANGES' C C Magic value C call dsa_get_flag_value('FLOAT',magic_float,status) if (status .ne. 0) go to 500 C C First initialise the output data array C totspace = nz*nkeep*dsa_typesize('float',status) call gen_fill(totspace,0,data) C C Initialise used spectra array C do i = 1,nkeep used(1,i) = 0 used(2,i) = 0 used(3,i) = 0 end do C C Decide how much around the edge to skip for the smoothing box C ipad1 = boxsize/2 if (mod(boxsize,2) .eq. 0) then ipad2 = ipad1 - 1 else ipad2 = ipad1 end if C C Sum all the data over the boxsize in the desired region C n = 0 do i = ipad1+1,ny-ipad2 do j = ipad1+1,nx-ipad2 C C Continue so long as this pixel hasn't be flagged for skipping C if (mon(1,j,i) .ne. 0) then n = n + 1 C C Zero a counter which will count the number of contributions C at each pixel. C do iz = 1,nz nsum(iz) = 0 end do C C Loop through the smoothing box. C do ii = i-ipad1,i+ipad2 do jj = j-ipad1,j+ipad2 C C Test each pixel in the spectral direction. If it's ok, then C add it into the data sum. If not then just count it. C nbad = 0 do k = 1,nz idata = input(k,jj,ii) if (idata .gt. magic_float) then data(k,n) = data(k,n) + idata nsum(k) = nsum(k) + 1 else nbad = nbad + 1 end if end do C C If this whole spectrum was flagged values, then flag C the spectrum as not to be used C if (nbad .eq. nz) mon(1,jj,ii) = 0 end do end do end if C C Normalise the data by the number of contributions to it. C if (mon(1,j,i) .ne. 0) then mon(1,j,i) = n used(1,n) = j used(2,n) = i do k = 1,nz if (nsum(k) .ne. 0) then data(k,n) = data(k,n)/float(nsum(k)) C C Create an error array depending on the error type C specified by the user C if (iwt .eq. 1) then d = data(k,n) err(k,n) = 1.0/sqrt(abs(d+1)) else if (iwt .eq. 2) then err(k,n) = 1.0 end if C C If a pixel has no contributions to it, then flag it C else data(k,n) = 0.0 err(k,n) = 100.0 end if end do end if end do end do 500 end subroutine fitt(nx,ny,nz,nkeep,mon,used,data,nco,ncomax,fit_type, : how_guess,fsr,alam0,dlam,niter,picture,params,err,lamda,covar, : alpha,lista,coefs,errcoefs) C+-------------------------------------------------------------------------------------- C FIT --- Routine to supervise the actual fitting of the lines. C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) NX: (integer) X dimension of the data C (>) NY: (integer) Y dimension of the data C (>) NZ: (integer) Z dimension of the data C (>) NKEEP: (integer) Total number of spectra to be used C (>) MON: (integer array MON(4,NX,NY)) Input 'monitor' array C (>) USED: (integer array USED(3,NKEEP)) Flags for used spectra C (>) DATA: (real array DATA(NZ,NKEEP)) Input data array C (>) NCO: (integer) Number of coefficients in the fit C (>) NCOMAX: (integer) Maximum number of coefficients C (>) FIT_TYPE: (integer) Flag for the fitting function to be used C (>) HOW_GUESS: (integer) Flag for the guessing algorithm to be used C (>) FSR: (real) Free spectral range of data C (>) ALAM0: (real) Wavelength zeropoint C (>) DLAM: (real) Wavelength increment C (>) NITER: (integer) Maximum number of iterations C (!) PICTURE: (logical) TRUE if a 2-d plot is currently on the screen C (>) PARAMS: (real array PARAMS(NCOMAX+4,NX,NY)) Final output array C (>) ERR: (real array ERR(NZ,NKEEP)) Data error array C (W) LAMDA: (real array LAMDA(NZ)) Workspace for spectral axis info C (W) COVAR: (real array COVAR(NCO,NCO)) Workspace for covariance matrix C (W) ALPHA: (real array ALPHA(NCO,NCO)) Workspace for curvarture matrix C (W) LISTA: (integer array LISTA(NCO)) Flags which parameters to be changed C (<) COEFS: (real array COEFS(NCO,NKEEP)) Output Coefs C (<) ERRCOEFS: (real array ERRCOEFS(NCO,NKEEP)) Error in output coefs C C C AUTHOR: Jim Lewis -- RGO C DATE: 14-SEP-1990 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nx,ny,nz,nkeep,nco,niter,how_guess,ncomax,fit_type integer lista(nco),mon(4,nx,ny),used(3,nkeep) real data(nz,nkeep),coefs(nco,nkeep),err(nz,nkeep),dlam, : covar(nco,nco),errcoefs(nco,nkeep),alpha(nco,nco),alam0, : lamda(nz),fsr,params(ncomax+4,nx,ny) logical picture C C Functions C integer dyn_element integer dsa_typesize integer ich_len C C Local variables C integer nfail,i,j,bytes_float,bytes_int,totspace,address, : slot,iwpt1,iwpt2,iwpt3,iwpt4,iwpt5,iwpt6,status,nspec,nrun, : np,wptr,tmon,k,slot2,l,nfail0,nrun0,iwpt7 real zmin,zmax,xvstart,yvstart,margin,deltax,deltay,xdev,ydev character message*80 C C Common block for dynamic memory. C include 'DYNAMIC_MEMORY' C C Numeric ranges C include 'NDP_NUMERIC_RANGES' C C A few initialisations for the fitting routine... C do i = 1,nco lista(i) = i end do nfail = 0 nrun = 0 status = 0 C C Get workspace for fitting routines... C bytes_float = dsa_typesize('float',status) bytes_int = dsa_typesize('int',status) totspace = 4*nco*bytes_float + 3*nco*bytes_int : + 2*nz*bytes_float call dsa_get_workspace(totspace,address,slot,status) if (status .ne. 0) go to 500 iwpt1 = dyn_element(address) iwpt2 = iwpt1 + nco*bytes_float iwpt3 = iwpt2 + nco*bytes_float iwpt4 = iwpt3 + nco*bytes_float iwpt5 = iwpt4 + nco*bytes_float iwpt6 = iwpt5 + 3*nco*bytes_int iwpt7 = iwpt6 + nz*bytes_float C C If guessing then start by setting up screen... C if (how_guess .ge. 3) then call clear_screen(status) if (status .ne. 0) then status = 0 go to 500 end if picture = .false. call pre_plot(nz,lamda,zmin,zmax,0.0,1.0,nkeep,xvstart,yvstart, : margin,deltax,deltay,xdev,ydev,np,slot2,wptr,status) if (status .ne. 0) then status = 0 go to 500 end if if (slot2 .ne. 0) call dsa_free_workspace(slot2,status) end if tmon = 0 C C Now do the fit for each spectrum C do l = 1,nkeep i = used(1,l) j = used(2,l) if (mon(1,i,j) .ne. 0) then nspec = mon(1,i,j) status = 0 C C Fill in the Z axis C do k = 1,nz lamda(k) = float(k) end do C C Take a guess at the coefs C call guess(how_guess,nkeep,nspec,data(1,nspec),used,lamda, : nx,ny,nz,fsr,nco,ncomax,params,alam0,dlam,tmon,xvstart, : yvstart,np,deltax,deltay,margin,zmin,zmax,coefs(1,nspec), : status) C C Check to see if the line is near an edge C call fiddle(coefs(1,nspec),data(1,nspec),err(1,nspec), : nz,fsr,dynamic_mem(iwpt6),dynamic_mem(iwpt7),lamda,status) C C Now do the fit C nfail0 = nfail nrun0 = nrun call sub_fit(niter,nz,lamda,dynamic_mem(iwpt6), : dynamic_mem(iwpt7),nco,lista,covar,alpha,fit_type,alam0, : dlam,dynamic_mem(iwpt1),dynamic_mem(iwpt2), : dynamic_mem(iwpt3),dynamic_mem(iwpt4),dynamic_mem(iwpt5), : coefs(1,nspec),errcoefs(1,nspec),nfail,nrun,status) end if if (nfail .gt. nfail0) then used(3,l) = -3 else used(3,l) = 1 end if if (mod(l,100) .eq. 0) then write(message,100) l,nkeep 100 format(i7,' spectra done out of',i7) call dsa_wruser(message(:ich_len(message))) call dsa_wrflush end if end do C C Send a message to the user... C write(message,101) nfail,niter 101 format(i8,' spectra failed to converge within',i3,' iterations') call dsa_wruser(message(:ich_len(message))) call dsa_wrflush write(message,102) nrun 102 format(i8,' spectra had runaway solutions') call dsa_wruser(message(:ich_len(message))) call dsa_wrflush call dsa_free_workspace(slot,status) 500 end subroutine sub_fit(niter,nz,lamda,data,err,nco,lista,covar, : alpha,fit_type,alam0,dlam,work1,work2,work3,work4,work5,coefs, : errcoefs,nfail,nrun,status) C+-------------------------------------------------------------------------------------- C SUB_FIT --- Routine to do the fitting of a single spectrum C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) NITER: (integer) Maximum number of iterations C (>) NZ: (integer) Z dimension of the data C (>) LAMDA: (real array LAMDA(NZ)) Spectral axis info C (>) DATA: (real array DATA(NZ)) Input data array C (>) ERR: (real array ERR(NZ)) Data error array C (>) NCO: (integer) Number of coefficients in the fit C (W) LISTA: (integer array LISTA(NCO)) Flags which parameters to be changed C (W) COVAR: (real array COVAR(NCO,NCO)) Workspace for covariance matrix C (W) ALPHA: (real array ALPHA(NCO,NCO)) Workspace for curvarture matrix C (>) FIT_TYPE: (integer) Flag for the fitting function to be used C (>) ALAM0: (real) Wavelength zeropoint C (>) DLAM: (real) Wavelength increment C (W) WORK1: (real array WORK1(NCO)) Workspace for fitting routine. C (W) WORK2: (real array WORK2(NCO)) Workspace for fitting routine. C (W) WORK3: (real array WORK3(NCO)) Workspace for fitting routine. C (W) WORK4: (real array WORK4(NCO)) Workspace for fitting routine. C (W) WORK5: (integer array WORK5(3*NCO)) Workspace for fitting routine. C (<) COEFS: (real array COEFS(NCO)) Output Coefs C (<) ERRCOEFS: (real array ERRCOEFS(NCO)) Error in output coefs C (<) NFAIL: (integer) Number of failed fits. C (<) NRUN: (integer) Number of runaway fits. C (<) STATUS: (integer) Status flag C C C AUTHOR: Jim Lewis -- RGO C DATE: 14-SEP-1990 C UPDATE: 21-DEC-1990 -- Modified error estimates to scale to RMS. (JRL) C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer niter,nz,nco,nfail,status,nrun,fit_type integer lista(nco),work5(3*nco) real lamda(nz),data(nz),err(nz),covar(nco,nco),alpha(nco,nco), : alam0,dlam,work1(nco),work2(nco),work3(nco),work4(nco), : coefs(nco),errcoefs(nco) C C External routines C external gauss,cauchy C C Local variables C integer ntriv,iter,ilam,j,ntriv_min,fit_status,status real alamda,olamda,chisq,dchi_max,sumwts,rms,magic_float logical again,run C C Common block with convergence parameters C common /conv/ dchi_max,ntriv_min C C Inherited status C if (status .ne. 0) then nfail = nfail + 1 go to 500 end if C C Initialise a few things... C iter = 0 alamda = -1.0 olamda = alamda again = .true. ntriv = 0 fit_status = 0 run = .false. call dsa_get_flag_value('FLOAT',magic_float,status) if (status .ne. 0) go to 500 C C Sum up the weights C sumwts = 0.0 do j = 1,nz sumwts = sumwts + 1.0/(err(j))**2 end do C C Iterate the fit until it converges or you run out of iterations... C do while ((iter .lt. niter) .and. (again) .and. : (fit_status .eq. 0)) iter = iter + 1 if (fit_type .eq. 1) then call mrqmin(lamda,1,data,err,nz,coefs,nco,lista,nco,covar, : alpha,nco,chisq,gauss,alamda,work1,work2,work3,work4,work5, : fit_status) else if (fit_type .eq. 2) then call mrqmin(lamda,1,data,err,nz,coefs,nco,lista,nco,covar, : alpha,nco,chisq,cauchy,alamda,work1,work2,work3,work4,work5, : fit_status) end if ilam = nint(olamda/alamda) call tau_converge(iter,ilam,alamda,chisq,ntriv,ntriv_min, : dchi_max,run,again) olamda = alamda end do C C Check to see whether the solution ever converged...If it didn't then C put magic values in for the coeffs and their errors C if (((iter .eq. niter) .and. (again)) .or. : (fit_status .ne. 0)) then nfail = nfail + 1 do j = 1,nco coefs(j) = magic_float errcoefs(j) = magic_float end do status = 1 C C If it did converge, then call the fitting routine once more to determine C the errors in the coeffs. C else alamda = 0.0 if (fit_type .eq. 1) then call mrqmin(lamda,1,data,err,nz,coefs,nco,lista,nco,covar, : alpha,nco,chisq,gauss,alamda,work1,work2,work3,work4,work5, : fit_status) else if (fit_type .eq. 2) then call mrqmin(lamda,1,data,err,nz,coefs,nco,lista,nco,covar, : alpha,nco,chisq,cauchy,alamda,work1,work2,work3,work4,work5, : fit_status) end if if (fit_status .ne. 0) status = 1 rms = sqrt(chisq/sumwts) do j = 1,nco if (covar(j,j) .lt. 0.0) then errcoefs(j) = -rms*sqrt(abs(covar(j,j))) else errcoefs(j) = rms*sqrt(covar(j,j)) end if end do C C Scale the results into wavelength units... C coefs(2) = alam0 + dlam*coefs(2) coefs(3) = abs(dlam)*coefs(3) errcoefs(2) = abs(dlam)*errcoefs(2) errcoefs(3) = abs(dlam)*errcoefs(3) C C If this was a runaway solution, then count it... C if (run) nrun = nrun + 1 end if 500 end subroutine limits(xmin,xmax,ymin,ymax,nx,ny,mask,thresh,xaxis, : yaxis,picture,mptr) C+-------------------------------------------------------------------------------------- C LIMITS --- Routine to set the limits on the domains the user C wants to restrict C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (<) XMIN: (real) Minimum value of X for consideration C (<) XMAX: (real) Maximum value of X for consideration C (<) YMIN: (real) Minimum value of Y for consideration C (<) YMAX: (real) Maximum value of Y for consideration C (>) NX: (integer) X dimension of the data C (>) NY: (integer) Y dimension of the data C (>) MASK: (logical) Instruction to use a mask image C (<) THRESH: (real array THRESH(2)) Min and max mask thresholds C (>) XAXIS: (real array XAXIS(NX)) X axis data array C (>) YAXIS: (real array YAXIS(NY)) Y axis data array C (>) PICTURE: (logical) Instruction to do a 2-d plot C (<) MPTR: (integer) Dynamic memory pointer for mask image. C C C AUTHOR: Jim Lewis -- RGO C DATE: 14-SEP-1990 C UPDATE: 19-DEC-1990 -- Took out restriction in spectral direction (JRL) C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nx,ny,mptr real xmin,xmax,ymin,ymax,thresh(2),xaxis(nx),yaxis(ny) logical mask,picture C C Functions C integer dyn_element logical par_batch logical par_quest C C Local variables C integer ndim,dims(2),nelm,status,slot,address real vmin(2),vmax(2),start(2),end(2) logical mapped,badpix,use_curs,batch character type_mask*8 data mapped /.false./ C C Numeric ranges C include 'NDP_NUMERIC_RANGES' C C Common block for dynamic memory. C include 'DYNAMIC_MEMORY' C C If using a mask, then open it and map it for later use C batch = par_batch() if (mask) then if (.not. mapped) then if (.not. batch) call par_cnpar('name_mask') call dsa_input('name_mask','name_mask',status) call ndp_get_image_info('name_mask',.true.,.false., : type_mask,badpix,status) call dsa_data_size('name_mask',2,ndim,dims,nelm,status) if (status .ne. 0) go to 500 C C If the mask image is not the same size as the input cube in C the x,y plane then it isn't going to be much use. C if (ndim .ne. 2) then call dsa_wruser('Mask image must be 2-d') call dsa_wrflush call dsa_close_structure('name_mask',status) go to 500 end if if ((dims(1) .ne. nx) .and. (dims(2) .ne. ny)) then call dsa_wruser : ('Mask image must be same size as data cube') call dsa_wrflush call dsa_close_structure('name_mask',status) go to 500 end if if (badpix) call dsa_use_flagged_values('name_mask',status) call dsa_map_data('name_mask','read','float',address,slot, : status) mptr = dyn_element(address) mapped = .true. end if C C Get the mask thresholds to be used. C if (.not. batch) call par_cnpar('mask_thresh') call par_rdary('mask_thresh',min_float,max_float,'a',' ', : 2,2,thresh) end if C C If you're using a 2-d plot, try and get the limits from the cursor C if (picture) then use_curs = par_quest('Use cursor to identify region?', : .true.) if (use_curs) call curse_limits(nx,ny,xaxis,yaxis,xmin,xmax, : ymin,ymax,status) end if C C If you're not using a 2-d plot, or if you are and not using the C cursor to define your region, or if the 2-d plot failed for any C reason, then input the limits by hand. C if ((.not. picture) .or. ((picture) .and. (.not. use_curs)).or. : (status .ne. 0)) then vmin(1) = 1.0 vmin(2) = 1.0 vmax(1) = float(nx) vmax(2) = float(ny) if (.not. batch) then call par_cnpar('start') call par_cnpar('end') end if call ndp_par_rdary('start',vmin,vmax,' ','pixels',2,2,start) call ndp_par_rdary('end',vmin,vmax,' ','pixels',2,2,end) xmin = start(1) ymin = start(2) xmax = end(1) ymax = end(2) if (status .ne. 0) status = 0 end if C C Draw a box of the defined limits C if (picture) then call pgwindow(1.0,float(nx),1.0,float(ny)) call pgmove(xmin,ymin) call pgdraw(xmax,ymin) call pgdraw(xmax,ymax) call pgdraw(xmin,ymax) call pgdraw(xmin,ymin) end if 500 end subroutine output(coefs,errcoefs,nx,ny,nco,ncomax,nkeep,mon, : used,fit_type,boxsize,params,errpar) C+-------------------------------------------------------------------------------------- C OUTPUT --- Routine to transfer the coefs to the output map C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) COEFS: (real array COEFS(NCO,NKEEP)) Input coefs C (>) ERRCOEFS: (real array ERRCOEFS(NCO,NKEEP)) Error in input coefs C (>) NX: (integer) X dimension of the data C (>) NY: (integer) Y dimension of the data C (>) NCO: (integer) Number of coefficients in the fit C (>) NCOMAX: (integer) Maximum number of coefficients in the fit C (>) NKEEP: (integer) Total number of spectra to be used C (>) MON: (integer array MON(4,NX,NY)) Input 'monitor' array C (>) USED: (integer array USED(3,NKEEP)) Flags for used spectra C (>) FIT_TYPE: (integer) Flag for fitting function. C (>) BOXSIZE: (integer) Size of binning box C (<) PARAMS: (real array PARAMS(NCOMAX+4,NX,NY)) Output data array C (<) ERRPAR: (real array ERRPAR(NCOMAX+4,NX,NY)) Output error array C C C AUTHOR: Jim Lewis -- RGO C DATE: 14-SEP-1990 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nx,ny,nkeep,nco,ncomax,fit_type,boxsize integer mon(4,nx,ny),used(3,nkeep) real coefs(nco,nkeep),errcoefs(nco,nkeep),params(ncomax+4,nx,ny), : errpar(ncomax+4,nx,ny) C C Local variables C integer i,j,k,l,nspec C C Transfer the results over... C do l = 1,nkeep i = used(1,l) j = used(2,l) if (mon(1,i,j) .gt. 0) then nspec = mon(1,i,j) do k = 1,nco params(k,i,j) = coefs(k,nspec) errpar(k,i,j) = errcoefs(k,nspec) end do mon(2,i,j) = fit_type mon(3,i,j) = used(3,l) mon(4,i,j) = boxsize end if end do end subroutine gauss(arg,nx,par,val,dval_dpar,npar) C+-------------------------------------------------------------------------------------- C GAUSS --- Routine to supply fitting function and Jacobean for C gaussian. Required for Press et. al. fitting routines. C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) ARG: (real array ARG(NX)) Place where function is evaluated C (>) NX: (integer) Number of independent variables C (>) PAR: (real array PAR(NPAR)) Most recent guess at fitting parameters C (<) VAL: (real) Value of function at ARG. C (<) DVAL_DPAR: (real array DVAL_DPAR(NPAR)) Value of Jacobean at ARG C (>) NPAR: (integer) Number of parameters in fit. C C C AUTHOR: Jim Lewis -- RGO C DATE: 14-SEP-1990 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer npar,nx real arg(nx),par(npar),val,dval_dpar(npar) C C Local variables C integer m real exparg,dx,sig2,expon,fsr,xx C C Common block with FSR C common /free/ fsr C C Evaluate the function and its Jacobean C sig2 = par(3)*par(3) xx = arg(1) m = nint((xx - par(2))/fsr) dx = xx - par(2) - float(m)*fsr exparg = dx*dx*0.5/sig2 if (exparg .ge. 8.0) then expon = 0.0 else expon = exp(-exparg) end if val = par(1)*expon + par(4) ! par(1) = A, par(2) = Xo, dval_dpar(1) = expon ! par(3) = sig, par(4) = back dval_dpar(2) = par(1)*dx*expon/sig2 dval_dpar(3) = par(1)*dx*dx*expon/(sig2*par(3)) dval_dpar(4) = 1.0 end subroutine cauchy(arg,nx,par,val,dval_dpar,npar) C+-------------------------------------------------------------------------------------- C CAUCHY --- Routine to supply fitting function and Jacobean for C Cauchy function. Required for Press et. al. fitting routines. C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) ARG: (real array ARG(NX)) Place where function is evaluated C (>) NX: (integer) Number of independent variables C (>) PAR: (real array PAR(NPAR)) Most recent guess at fitting parameters C (<) VAL: (real) Value of function at ARG. C (<) DVAL_DPAR: (real array DVAL_DPAR(NPAR)) Value of Jacobean at ARG C (>) NPAR: (integer) Number of parameters in fit. C C C AUTHOR: Jim Lewis -- RGO C DATE: 30-NOV-1990 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer npar,nx real arg(nx),par(npar),val,dval_dpar(npar) C C Local variables C integer m real fsr,denom,sig2,xx,dx,fac C C Common block with FSR C common /free/ fsr C C Evaluate the function and its Jacobean C sig2 = par(3)*par(3) xx = arg(1) m = nint((xx - par(2))/fsr) dx = xx - par(2) - float(m)*fsr fac = 0.5*dx*dx/sig2 denom = 1.0/(1.0 + fac) val = par(1)*denom + par(4) dval_dpar(1) = denom dval_dpar(2) = par(1)*dx*denom*denom/sig2 dval_dpar(3) = par(1)*2.0*fac*denom*denom/par(3) dval_dpar(4) = 1.0 end subroutine guess(how_guess,nkeep,nspec,data,used,lamda,nx,ny, : nz,fsr,nco,ncomax,params,alam0,dlam,tmon,xvstart,yvstart, : np,deltax,deltay,margin,zmin,zmax,coefs,status) C+-------------------------------------------------------------------------------------- C GUESS --- Routine to supervise the guessing algorithms. C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) HOW_GUESS: (integer) Flag for the guessing algorithm to be used C (>) NKEEP: (integer) Total number of spectra to be used C (>) NSPEC: (integer) Which spectrum in input data to use. C (>) DATA: (real array DATA(NZ)) Input data array C (>) USED: (integer array USED(3,NKEEP)) Flags for used spectra C (W) LAMDA: (real array LAMDA(NZ)) Workspace for spectral axis info C (>) NX: (integer) X dimension of the data C (>) NY: (integer) Y dimension of the data C (>) NZ: (integer) Z dimension of the data C (>) FSR: (real) Free spectral range of data C (>) NCO: (integer) Number of coefficients in the fit C (>) NCOMAX: (integer) Maximum number of coefficients C (>) PARAMS: (real array PARAMS(NCOMAX+4,NX,NY)) Final output array C (>) ALAM0: (real) Wavelength zeropoint C (>) DLAM: (real) Wavelength increment C (!) TMON: (integer) How many plots are on the page C (>) XVSTART: (real) Where in the viewport to put the next plot C (>) YVSTART: (real) Where in the viewport to put the next plot C (>) NP: (integer) Number of plots on a side C (>) DELTAX: (real) X size of each plot C (>) DELTAY: (real) Y size of each plot C (>) MARGIN: (real) Size of the margin for each plot C (>) ZMIN: (real) Minimum data value C (>) ZMAX: (real) Maximum data value C (<) COEFS: (real array COEFS(NCO)) Output guessed coefs C (!) STATUS: (integer) Status variable C C C AUTHOR: Jim Lewis -- RGO C DATE: 14-SEP-1990 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nkeep,nz,nco,ncomax,how_guess,nspec,nx,ny,tmon,np,status integer used(3,nkeep) real data(nz),lamda(nz),coefs(nco),fsr, : params(ncomax+4,nx,ny),alam0,dlam,xvstart,yvstart,deltax,deltay, : margin,zmin,zmax C C Functions C integer dsa_typesize integer dyn_element C C Local variables C integer totspace,bytes_float,slot,wpt1,wpt2, : address,nrem,ndiv,iadd,guess_status real xv1,yv1,xv2,yv2,dmin,dmax logical quick C C Common block for dynamic memory. C include 'DYNAMIC_MEMORY' C C Guess by looking at moments C if (how_guess .eq. 1) then call guess_moms(data,nz,fsr,coefs) C C Use a previous result as guess C else if (how_guess .eq. 2) then call guess_previous(data,nx,ny,nz,fsr,nco,ncomax, : nspec,nkeep,used,params,alam0,dlam,coefs) C C Do the guess by hand. C else if ((how_guess .eq. 3) .or. (how_guess .eq. 4)) then C C Find where on the viewport to put this plot C tmon = tmon + 1 if (np .eq. 1) then xv1 = xvstart yv1 = yvstart else nrem = mod(tmon,np) ndiv = (tmon - 1)/np if (nrem .eq. 1) then xv1 = xvstart yv1 = yvstart + float(ndiv)*(deltay + 2*margin) iadd = 0 else iadd = iadd + 1 xv1 = xvstart + float(iadd)*(deltax + 2*margin) end if end if xv2 = xv1 + deltax yv2 = yv1 + deltay C C Get some workspace C bytes_float = dsa_typesize('float',status) totspace = 2*nz*bytes_float call dsa_get_workspace(totspace,address,slot,status) if (status .ne. 0) go to 500 wpt1 = dyn_element(address) wpt2 = wpt1 + nz*bytes_float C C Now do the plot C if ((np .eq. 1) .or. (mod(tmon,np*np) .eq. 1)) call pgadvance call pgvsize(xv1,xv2,yv1,yv2) call minmax(nz,data,dmin,dmax) call pgwindow(zmin,zmax,dmin,dmax) call pgbox('BCNT',0.0,0,'BCNTS',0.0,0) call ndp_pgbin(lamda,data,dynamic_mem(wpt1), : dynamic_mem(wpt2),nz,.false.,.true.) C C Do the guessing... C quick = .true. if (how_guess .eq. 3) quick = .false. call guess_by_hand(quick,data,nz,lamda,dynamic_mem(wpt1),coefs, : guess_status) call dsa_free_workspace(slot,status) if (guess_status .ne. 0) status = 1 end if 500 end subroutine fiddle(coefs_guess,datain,err,nz,fsr,datout,errout, : lamda,status) C+-------------------------------------------------------------------------------------- C FIDDLE --- Routine to fiddle the input data array if the line is C too close to either end of the spectrum C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) COEFS_GUESS: (real array COEFS_GUESS(4)) The best guess at the fit parameters C (>) DATAIN: (real array DATAIN(NZ)) Input data array C (>) ERR: (real array ERR(NZ,NKEEP)) Data error array C (>) NZ: (integer) Z dimension of the data C (>) FSR: (real) Free spectral range of data C (<) DATOUT: (real array DATOUT(NZ)) Output data array C (<) ERROUT: (real array ERROUT(NZ)) Output error array C (<) LAMDA: (real array LAMDA(NZ)) Workspace for spectral axis info C (!) STATUS: (integer) Status variable C C C AUTHOR: Jim Lewis -- RGO C DATE: 14-SEP-1990 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nz,status real coefs_guess(4),datain(nz),fsr,datout(nz),lamda(nz),err(nz), : errout(nz) C C Local variables C integer nfsr,ist,isub,index,i,nfsr2 real di1,di2,div1,div2,fsr2 C C Inherited status C if (status .ne. 0) go to 500 C C First find out where the line is w.r.t. both edges. C di1 = coefs_guess(2) di2 = float(nz) - coefs_guess(2) + 1.0 C C Now find out if it is too close to an edge (defined as being C within 1 sigma...) C div1 = di1/coefs_guess(3) div2 = di2/coefs_guess(3) if ((div1 .ge. 1.0) .and. (div2 .ge. 1.0)) then isub = 0 else if ((div1 .lt. 1.0) .and. (div2 .lt. 1.0)) then isub = 0 else if ((div1 .ge. 1.0) .and. (div2 .lt. 1.0)) then isub = -1 else if ((div1 .lt. 1.0) .and. (div2 .ge. 1.0)) then isub = 1 end if C C If the line is well centred then just create a wavelength C array and a copy of the old data (isub = 0). If the line C runs off either side, then pad the line with half an FSR's C worth of data from a region 1 FSR away. Set up wavelengths C appropriately. C fsr2 = 0.5*fsr nfsr2 = nint(fsr2) nfsr = nint(fsr) do i = 1,nz ist = i - nfsr2*isub if (ist .lt. 1) then index = ist + nfsr else if (ist .gt. nz) then index = ist - nfsr else index = ist end if lamda(i) = float(ist) datout(i) = datain(index) errout(i) = err(index) end do 500 end subroutine guess_moms(data,nz,fsr,coefs) C+-------------------------------------------------------------------------------------- C GUESS_MOMS --- Guess at the parameters by looking at the moments in the spectrum C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) DATA: (real array DATA(NZ)) Input data array C (>) NZ: (integer) Z dimension of the data C (>) FSR: (real) Free spectral range C (<) COEFS: (real array COEFS(4)) The best guess at the fit parameters C C C AUTHOR: Jim Lewis -- RGO C DATE: 14-SEP-1990 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nz real data(nz),coefs(4),fsr C C Local variables C integer i,irem,i1,i2,di,di2,ist,ifn,ifsr,itest,ilook real dmin,dmax,sum,sumsq,l,d,ave,sig,tot,sig2 C C Numeric Ranges C include 'NDP_NUMERIC_RANGES' C C Find the peak in the data C dmax = min_float dmin = max_float do i = 1,nz dmin = min(dmin,data(i)) if (data(i) .gt. dmax) then dmax = data(i) irem = i end if end do ifsr = nint(fsr) C C Subtract 'background' C dmax = dmax - dmin C C Find half max points C i1 = irem itest = i1 ilook = 0 do while ((data(itest)-dmin .gt. 0.5*dmax) .and. (ilook .lt. nz)) ilook = ilook + 1 i1 = i1 + 1 if (i1 .gt. nz) then itest = i1 - ifsr else itest = i1 end if end do i2 = irem itest = i2 ilook = 0 do while ((data(itest)-dmin .gt. 0.5*dmax) .and. (ilook .lt. nz)) ilook = ilook + 1 i2 = i2 - 1 if (i2 .lt. 1) then itest = i2 + ifsr else itest = i2 end if end do C C Define crude line width and box over which to sum data C di = (i1 - i2)/2 di2 = 2*di ist = irem - di2 ifn = irem + di2 C C Now find mean and sigma C sum = 0.0 sumsq = 0.0 tot = 0.0 do i = ist,ifn if (i .lt. 1) then itest = i + ifsr else if (i .gt. nz) then itest = i - ifsr else itest = i end if l = float(i) d = data(itest) - dmin sum = sum + l*d sumsq = sumsq + l*l*d tot = tot + d end do ave = sum/tot sig2 = (sumsq/tot - ave*ave) if (sig2 .lt. 0.0) then sig = abs(sig2) else sig = sqrt(sumsq/tot - ave*ave) end if C C Fill in the coefs C coefs(1) = dmax ! Peak coefs(2) = ave ! Position coefs(3) = sig ! Width coefs(4) = dmin ! Background end subroutine guess_previous(data,nx,ny,nz,fsr,nco,ncomax,nspec, : nkeep,used,params,alam0,dlam,coefs) C+-------------------------------------------------------------------------------------- C GUESS_PREVIOUS --- Guess at the parameters by looking at a previous fit C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) DATA: (real array DATA(NZ)) Input data array C (>) NX: (integer) X dimension of the data C (>) NY: (integer) Y dimension of the data C (>) NZ: (integer) Z dimension of the data C (>) FSR: (real) Free spectral range C (>) NCO: (integer) Number of coefs C (>) NCOMAX: (integer) Maximum number of coefs C (>) NKEEP: (integer) Number of spectra kept C (>) NSPEC: (integer) Number of the spectrum C (>) USED: (integer array USED(3,NKEEP)) Flags for used spectra C (>) PARAMS: (real array PARAMS(NCOMAX+4,NX,NY)) Final output array C (>) ALAM0: (real) Wavelength zeropoint C (>) DLAM: (real) Wavelength increment C (<) COEFS: (real array COEFS(NCO)) Best guess at coefs C C C AUTHOR: Jim Lewis -- RGO C DATE: 14-SEP-1990 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nz,nx,ny,ncomax,nco,nspec,nkeep integer used(3,nkeep) real fsr,data(nz),coefs(nco),params(ncomax+4,nx,ny),alam0,dlam C C Local variables C integer i,j,k,status real magic_float C C Magic value C call dsa_get_flag_value('FLOAT',magic_float,status) if (status .ne. 0) go to 500 C C Load the old values C i = used(1,nspec) j = used(2,nspec) do k = 1,nco coefs(k) = params(k,i,j) end do C C If the old fit doesn't exist for one reason or another, then do C a moments guess... C if (coefs(1) .le. magic_float) then call guess_moms(data,nz,fsr,coefs) else coefs(1) = coefs(1) coefs(2) = (coefs(2) - alam0)/dlam coefs(3) = coefs(3)/abs(dlam) coefs(4) = coefs(4) end if 500 end subroutine guess_by_hand(quick,data,nz,lamda,work1,coefs, : status) C+-------------------------------------------------------------------------------------- C GUESS_BY_HAND --- Guess at the parameters by showing the program the salient C features of the spectrum with a cursor C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) QUICK: (logical) TRUE if the user is doing the quicker version C (>) DATA: (real array DATA(NZ)) Input data array C (>) NZ: (integer) Z dimension of the data C (>) LAMDA: (real array LAMDA(NZ)) Spectral axis info C (W) WORK1: (real array WORK1(NZ)) Workspace for plotting routines C (<) COEFS: (real array COEFS(4)) The best guess at the fit parameters C (!) STATUS: (integer) Status flag C C C AUTHOR: Jim Lewis -- RGO C DATE: 14-SEP-1990 C UPDATE: 09-JAN-1990 -- Added quick version. (JRL) C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nz,status real data(nz),coefs(4),lamda(nz),work1(nz) logical quick C C Functions C integer ich_len C C Maximum number of contiuum points allowed C integer max_cont parameter (max_cont=10) C C Local variables C integer ncont,nline,npeak,i,ixc,ist,ifn,irem real zmin,zmax,dmin,dmax,xcont(max_cont),xline(2),xpeak,sum, : halfpeak,slope,aint,halfwidth,back,fsr,magic_float character message*80 C C Common block with FSR C common /free/ fsr C C Redefine the window. C call dsa_get_flag_value('FLOAT',magic_float,status) if (status .ne. 0) go to 500 zmin = 1.0 zmax = float(nz) do i = 1,nz lamda(i) = float(i) end do call minmax(nz,data,dmin,dmax) call pgwindow(zmin,zmax,dmin,dmax) C C Choose the continuum points C if (.not. quick) then write(message,100) max_cont 100 format('You may choose up to',i4,' contiuum points') call dsa_wruser(' ') call dsa_wrflush call dsa_wruser(message(:ich_len(message))) call dsa_wrflush call tcursor(max_cont,zmin,zmax,dmin,dmax,nz,data,ncont, : xcont) if (ncont .eq. 0) then call dsa_wruser : ('You MUST choose at least one contiuum point!!!') call dsa_wrflush coefs(1) = magic_float go to 500 end if end if C C Mark out the limits of the line C call dsa_wruser(' ') call dsa_wrflush call dsa_wruser('Delimit the line with 1 point on either side') call dsa_wrflush call tcursor(2,zmin,zmax,dmin,dmax,nz,data,nline,xline) if (nline .lt. 2) then call dsa_wruser : ('You MUST choose 2 line delimiting points!!!') call dsa_wrflush coefs(1) = magic_float go to 500 end if C C If doing the quick version, then the line delimiters also delimit C the continuum C if (quick) then xcont(1) = xline(1) xcont(2) = xline(2) ncont = 2 end if C C Mark the peak of the line C call dsa_wruser(' ') call dsa_wrflush call dsa_wruser('Now show me where the peak is...') call dsa_wrflush call tcursor(1,zmin,zmax,dmin,dmax,nz,data,npeak,xpeak) if (npeak .eq. 0) then call dsa_wruser('You MUST show me where the peak is!!!') call dsa_wrflush coefs(1) = magic_float go to 500 end if C C Now set up the guesses based on what you've been told. C Obviously the peak has been pointed out and the background C is the average of the continuum points. The centre is the C position of the peak... C coefs(2) = xpeak sum = 0.0 do i = 1,ncont ixc = nint(xcont(i)) sum = sum + data(ixc) end do back = sum/float(ncont) coefs(4) = back coefs(1) = data(nint(xpeak)) - back C C Now get an estimate of the width by finding where the data falls to C half the value of the peak... C halfpeak = 0.5*coefs(1) + back ist = min(nint(xline(1)),nint(xline(2))) ifn = max(nint(xline(1)),nint(xline(2))) do i = ist,nint(xpeak) if (data(i) .le. halfpeak) irem = i end do slope = (data(irem+1) - data(irem)) aint = data(irem) - slope*float(irem) sum = abs((halfpeak - aint)/slope - xpeak) do i = ifn,nint(xpeak),-1 if (data(i) .le. halfpeak) irem = i end do slope = -(data(irem-1) - data(irem)) aint = data(irem) - slope*float(irem) sum = sum + abs((halfpeak - aint)/slope - xpeak) halfwidth = sum*0.5 coefs(3) = halfwidth/sqrt(2.0*log(2.0)) C C Now plot up the resulting guess in red C call gauss_eval(coefs,4,lamda,nz,fsr,work1) call pgsci(2) call pgscr(2,1.0,0.0,0.0) call pgline(nz,lamda,work1) call pgsci(1) C C If the routine has failed somewhere, then fill in magic values... C 500 if (coefs(1) .eq. magic_float) then do i = 2,4 coefs(i) = magic_float end do status = 1 end if end subroutine guess_other(nx,ny,nz,np,monitor,xv1,yv1,deltax,deltay, : margin,zmin,zmax,data,nco,ncomax,nkeep,mon,used,params,alam0, : dlam,coefs,status) C+-------------------------------------------------------------------------------------- C GUESS_OTHER --- Get guess for one spectrum from the fit coefs of C another C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) NX: (integer) X dimension of the data C (>) NY: (integer) Y dimension of the data C (>) NZ: (integer) Z dimension of the data C (>) NP: (integer) Number of plots on a side C (>) MONITOR: (integer array MONITOR(3,NP*NP)) Plot monitor array C (>) XV1: (real) X position of lower left hand corner of screen C (>) YV1: (real) Y position of lower left hand corner of screen C (>) DELTAX: (real) X size of plot in inches C (>) DELTAY: (real) Y size of plot in inches C (>) MARGIN: (integer) Size of margin in inches C (>) ZMIN: (real) Minimum value of Z C (>) ZMAX: (real) Maximum value of Z C (>) DATA: (real array DATA(NZ)) Data array of last spectrum plotted C (>) NCO: (integer) Number of coefs in the fit C (>) NCOMAX: (integer) Maximum number of coefs allowed C (>) NKEEP: (integer) Total number of spectra to be used C (>) MON: (integer array MON(4,NX,NY)) Monitor array C (>) USED: (integer array USED(3,NKEEP)) Flags for used spectra C (>) PARAMS: (real array PARAMS(NCOMAX+4,NX,NY)) Final output array C (>) ALAM0: (real) Wavelength zeropoint C (>) DLAM: (real) Wavelength increment C (<) COEFS: (real array COEFS(NCO,NKEEP)) Output Coefs C (!) STATUS: (integer) Status variable C C C AUTHOR: Jim Lewis -- RGO C DATE: 10-JAN-1991 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nx,ny,nz,np,nco,ncomax,nkeep,status integer monitor(3,np*np),used(3,nkeep),mon(4,nx,ny) real xv1,yv1,deltax,deltay,margin,zmin,zmax,data(nz),coefs(nco), : params(ncomax+4,nx,ny),alam0,dlam C C Functions C integer ich_fold integer ich_len C C Local variables C integer ix,iy,npage,nspec,dumint real xx,yy logical ok character ch*1,message*80 C C Select the spectrum to use as template C status = 0 ok = .false. do while (.not. ok) call dsa_wruser : ('Choose spectrum to use as template (Q to quit)') call dsa_wrflush call pgcurse(xx,yy,ch) dumint = ich_fold(ch) if (ch .eq. 'Q') then status = 1 go to 500 end if C C Find out which cell on the page you're dealing with C call find_cell(np,nz,xx,yy,monitor,xv1,yv1,deltax,deltay, : margin,zmin,zmax,data,ix,iy,npage,nspec) if (mon(3,used(1,nspec),used(2,nspec)) .le. 0) then write(message,100) used(1,nspec),used(2,nspec) 100 format('Cell ',2i5, : ' doesn''t have a valid fit...try again...') call dsa_wruser(message(:ich_len(message))) call dsa_wrflush else ok = .true. end if end do C C Now copy over the fit coefs for that cell C coefs(1) = params(1,used(1,nspec),used(2,nspec)) coefs(2) = (params(2,used(1,nspec),used(2,nspec)) - alam0)/dlam coefs(3) = params(3,used(1,nspec),used(2,nspec))/abs(dlam) coefs(4) = params(4,used(1,nspec),used(2,nspec)) 500 end subroutine tcursor(max_pts,zmin,zmax,dmin,dmax,nz,data,npts,xpos) C+-------------------------------------------------------------------------------------- C CURSOR --- Use cursor to mark various things on a spectrum C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) MAX_PTS: (integer) Maximum number of points to be selected C (>) ZMIN: (real) Lower limit of spectral axis C (>) ZMAX: (real) Upper limit of spectral axis C (>) DMIN: (real) Lower limit of data axis C (>) DMAX: (real) Upper limit of data axis C (>) NZ: (integer) Z dimension of the data C (>) DATA: (real array DATA(NZ)) Input data array C (<) NPTS: (integer) Number of points eventually chosen C (<) XPOS: (real array xpos(MAX_PTS)) The X axis position of the chosen points C C C AUTHOR: Jim Lewis -- RGO C DATE: 14-SEP-1990 C UPDATE: 09-JAN-1991 -- Can now quit by putting cursor outside frame (JRL) C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer max_pts,npts,nz real zmin,zmax,dmin,dmax,xpos(max_pts),data(nz) C C Functions C integer ich_fold C C Local variables C integer idummy real xx,yy character ch C C Loop for the maximum number of points allowed or until the C character returned is 'Q'. Can also quit by putting cursor C outside of frame. C npts = 0 ch = ' ' call pgwindow(zmin,zmax,dmin,dmax) do while ((npts .lt. max_pts) .and. (ch .ne. 'Q')) call pgcurse(xx,yy,ch) idummy = ich_fold(ch) if ((xx .lt. zmin) .or. (xx .gt. zmax) .or. : (yy .lt. dmin) .or. (yy .gt. dmax)) ch = 'Q' if (ch .ne. 'Q') then npts = npts + 1 xpos(npts) = xx call pgpoint(1,xx,data(nint(xx)),3) end if end do end subroutine gauss_eval(coefs,nc,zaxis,nz,fsr,yaxis) C+-------------------------------------------------------------------------------------- C GAUSS_EVAL --- Evaluate the Gauss fit along the spectral axis for plots C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) COEFS: (real array COEFS(NC)) The fit parameters C (>) NC: (integer) Dimensions of coef array C (>) ZAXIS: (real array ZAXIS(NZ)) Spectral (X) axis values C (>) NZ: (integer) Z dimension of the data C (>) FSR: (real) Free spectral range C (<) YAXIS: (real array YAXIS(NZ)) Data (Y) axis values of fit C C C AUTHOR: Jim Lewis -- RGO C DATE: 14-SEP-1990 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nc,nz real coefs(nc),zaxis(nz),yaxis(nz),fsr C C Local variables C integer i,m real width,peak,back,cent,arg C C Set up some constants C width = 0.5/coefs(3)**2 peak = coefs(1) back = coefs(4) cent = coefs(2) C C Now evaluate all gaussians... C do i = 1,nz m = nint((zaxis(i)-cent)/fsr) arg = -width*(zaxis(i)-cent-float(m)*fsr)**2 if (abs(arg) .lt. 40.0) then yaxis(i) = peak*exp(arg) + back else yaxis(i) = back end if end do end subroutine cauchy_eval(coefs,nc,zaxis,nz,fsr,yaxis) C+-------------------------------------------------------------------------------------- C CAUCHY_EVAL --- Evaluate the Cauchy fit along the spectral axis for plots C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) COEFS: (real array COEFS(NC)) The fit parameters C (>) NC: (integer) Dimension of coefs array C (>) ZAXIS: (real array ZAXIS(NZ)) Spectral (X) axis values C (>) NZ: (integer) Z dimension of the data C (>) FSR: (real) Free spectral range of data C (<) YAXIS: (real array YAXIS(NZ)) Data (Y) axis values of fit C C C AUTHOR: Jim Lewis -- RGO C DATE: 30-NOV-1990 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nc,nz real coefs(nc),zaxis(nz),fsr,yaxis(nz) C C Local variables C integer i,m real peak,back,cent,width C C Set up some constants C peak = coefs(1) back = coefs(4) cent = coefs(2) width = 0.5/coefs(3)**2 C C Now evaluate the function C do i = 1,nz m = nint((zaxis(i)-cent)/fsr) yaxis(i) = peak/(1.0+width*(zaxis(i)-cent-float(m)*fsr)**2) : + back end do end subroutine minmax(n,data,dmin,dmax) C+-------------------------------------------------------------------------------------- C MINMAX --- Find the min and max of a data array C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) N: (integer) Number of elements in array C (>) DATA: (real array DATA(N)) Input data array C (<) DMIN: (real) Minimum value of array minus one C (<) DMAX: (real) Maximum value of array plus one C C C AUTHOR: Jim Lewis -- RGO C DATE: 14-SEP-1990 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer n real data(n),dmin,dmax C C Local variables C integer i real yrange(2) logical scale_y C C Numeric Ranges C include 'NDP_NUMERIC_RANGES' C C Common block with Y scaling info for line plotting routines C common /scaling/ scale_y,yrange C C If you have set ranges, then use them C if (scale_y) then dmin = yrange(1) dmax = yrange(2) C C Otherwise find min and max C else dmin = max_float dmax = min_float do i = 1,n dmin = min(dmin,data(i)) dmax = max(dmax,data(i)) end do C C Pad it out... C dmin = dmin - 0.5 dmax = dmax + 0.5 end if end subroutine clear_screen(status) C+-------------------------------------------------------------------------------------- C CLEAR_SCREEN --- Clear the plotting surface C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (!) STATUS: (integer) Status variable C C C AUTHOR: Jim Lewis -- RGO C DATE: 28-NOV-1990 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer status C C Local variables C integer length character pgstate*20,softdev*16 C C If the plot device is already open, then close it...then reopen C call pgqinf('state',pgstate,length) if (pgstate .eq. 'OPEN') call pgend call var_getchr('softdev',0,0,softdev,status) if (status .ne. 0) then call dsa_wruser : ('No screen device selected -- use SOFT command') call dsa_wrflush go to 500 end if call pgbegin(0,softdev,1,1) call pgask(.false.) 500 end subroutine two_d_plot(nx,ny,xaxis,yaxis,picture,status) C+-------------------------------------------------------------------------------------- C TWO_D_PLOT --- Do a colour 2-d plot C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) NX: (integer) X dimension of the data C (>) NY: (integer) Y dimension of the data C (>) XAXIS: (real array XAXIS(NX)) X axis data C (>) YAXIS: (real array YAXIS(NY)) Y axis data C (!) PICTURE: (logical) TRUE if a 2-d plot is on the screen C (!) STATUS: (integer) Status variable C C C AUTHOR: Jim Lewis -- RGO C DATE: 28-NOV-1990 C UPDATE: 28-JAN-1992 -- Modified to use new plotting routines (JRL) C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nx,ny,status real xaxis(nx),yaxis(ny) logical picture C C Functions C integer dyn_element integer dsa_typesize logical par_qnum logical par_quest C C Local variables C integer ndim,dims(10),nelm,address,slot,pptr,stapix(2), : endpix(2),ignore,wslot,wptr,bytes_int,totspace real start(2),end(2),low,high,oldlow,oldhigh character type*8,table*20,name_2d*80,oldname_2d*80,oldtable*20 logical badpix,open,same,alloc,init data open /.false./ data init /.true./ C C Numeric ranges C include 'NDP_NUMERIC_RANGES' C C Common block for dynamic memory. C include 'DYNAMIC_MEMORY' C C Open up... C if (init) then call clear_screen(status) init = .false. end if C C Do you want the same plot as before? C if (open) then same = par_quest('Display same image you had last time?', : .true.) if (.not. same) then if (badpix) : call dsa_post_process_flagged_values('name_2d',status) call dsa_unmap(slot,status) call dsa_close_structure('name_2d',status) if (status .ne. 0) go to 500 open = .false. if (alloc) then alloc = .false. call dsa_free_workspace(wslot,status) if (status .ne. 0) go to 500 end if end if endif C C If a file hasn't been opened yet, then get the file name. C if (.not. open) then call var_getchr('name_2d',0,0,oldname_2d,ignore) call par_qstr('Name of 2-d image',oldname_2d,.false.,.false., : name_2d) call var_setchr('name_2d',0,0,name_2d,ignore) call var_getchr('table',0,0,oldtable,ignore) call par_qstr('Name of lookup table',oldtable,.false.,.false., : table) call var_setchr('table',0,0,table,ignore) call dsa_named_input('name_2d',name_2d,status) C C Check and make sure the file sizes are the same as for the C data cube C call ndp_get_image_info('name_2d',.true.,.false., : type,badpix,status) call dsa_data_size('name_2d',10,ndim,dims,nelm,status) if (ndim .ne. 2) then call dsa_wruser('Plot image must be 2-d') call dsa_wrflush call dsa_close_structure('name_2d',status) go to 500 end if if ((dims(1) .ne. nx) .and. (dims(2) .ne. ny)) then call dsa_wruser : ('Plot image must be same size as data cube') call dsa_wrflush call dsa_close_structure('name_2d',status) go to 500 end if if (badpix) call dsa_use_flagged_values('name_2d',status) if (status .ne. 0) go to 500 C C Map the data, setup some variables for the plots and read the C look up table name C call dsa_map_data('name_2d','read','float',address,slot, : status) pptr = dyn_element(address) open = .true. C C Workspace for colour index array C bytes_int = dsa_typesize('int',status) totspace = bytes_int*nx*ny call dsa_get_workspace(totspace,address,wslot,status) if (status .ne. 0) go to 500 wptr = dyn_element(address) alloc = .true. end if C C Each time the plot is done, ask for min and max to plot C call var_getnum('low',0,0,oldlow,ignore) call var_getnum('high',0,0,oldhigh,ignore) ignore = par_qnum('Low value for plot',min_float,max_float, : oldlow,.true.,' ',low) ignore = par_qnum('High value for plot',min_float,max_float, : oldhigh,.true.,' ',high) if (low .eq. high) then call dsa_wruser('Plot min and max can''t be the same!') call dsa_wrflush go to 500 end if call var_setnum('low',0,0,low,ignore) call var_setnum('high',0,0,high,ignore) C C Create colour index array as it's needed C if ((.not. same) .or. (low .ne. oldlow) .or. : (high .ne. oldhigh)) then call ndp_image_index(nx*ny,low,high,dynamic_mem(pptr),badpix, : dynamic_mem(wptr),status) if (status .ne. 0) go to 500 end if C C Set up some axis info C stapix(1) = 1 stapix(2) = 1 endpix(1) = nx endpix(2) = ny start(1) = xaxis(stapix(1)) start(2) = yaxis(stapix(2)) end(1) = xaxis(endpix(1)) end(2) = yaxis(endpix(2)) C C Do the plot C call disp(dynamic_mem(wptr),nx,ny,stapix,endpix,start,end, : high,low,picture,table,status) 500 end subroutine disp(array,nx,ny,stapix,endpix,start,end,high, : low,picture,table,status) C+-------------------------------------------------------------------------------------- C DISP --- Do a standard colour 2-d plot C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) ARRAY (integer array ARRAY(NX,NY)) Colour index data C (>) NX: (integer) X dimension of the data C (>) NY: (integer) Y dimension of the data C (>) STAPIX: (integer array STAPIX(2)) Pixel positions of starting points C (>) ENDPIX: (integer array ENDPIX(2)) Pixel positions of ending points C (>) START: (real array START(2)) Axis values of starting points C (>) END: (real array END(2)) Axis values of ending points C (>) HIGH: (real) High value for plot C (>) LOW: (real) Low value for plot C (!) PICTURE: (logical) TRUE if a 2-d plot is on the screen C (>) TABLE: (character) LUT for plot C (!) STATUS: (integer) Status variable C C C AUTHOR: Jim Lewis -- RGO C DATE: 25-FEB-1991 C UPDATE: 29-JAN-1992 -- Modified to use new NDP_ plotting routines (JRL) C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nx,ny,stapix(2),endpix(2),status integer array(nx,ny) real start(2),end(2),high,low character table*(*) logical picture C C Local variables C real ximv(2),yimv(2),square C C Clear the screen and plot the array C call clear_screen(status) picture = .false. call ndp_image_lut(table,status) if (status .ne. 0) go to 500 call ndp_image_viewport(stapix,endpix,1.0,'L',ximv,yimv, : square) call ndp_image_plot(array,nx,ny,stapix,endpix, : start,end,high,low,' ','AR',ximv,yimv) picture = .true. 500 end subroutine curse_limits(nx,ny,xaxis,yaxis,xmin,xmax,ymin,ymax, : status) C+-------------------------------------------------------------------------------------- C CURSE_LIMITS --- Get limits by using the cursor C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) NX: (integer) X dimension of the data C (>) NY: (integer) Y dimension of the data C (>) XAXIS: (real array XAXIS(NX)) X axis data C (>) YAXIS: (real array YAXIS(NY)) Y axis data C (<) XMIN: (real) Minimum value of X for consideration C (<) XMAX: (real) Maximum value of X for consideration C (<) YMIN: (real) Minimum value of Y for consideration C (<) YMAX: (real) Maximum value of Y for consideration C (<) STATUS: (integer) Status variable C C C AUTHOR: Jim Lewis -- RGO C DATE: 28-NOV-1990 C UPDATE: 19-DEC-1990 -- Modified to call routine CURSOR_POINT. (JRL) C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nx,ny,status real xaxis(nx),yaxis(ny),xmin,xmax,ymin,ymax C C Local variables C real xid(2),yid(2) logical quit1,quit2 C C Choose the points... C call cursor_point(nx,ny,xaxis,yaxis,xid(1),yid(1),quit1) if (.not. quit1) : call cursor_point(nx,ny,xaxis,yaxis,xid(2),yid(2),quit2) C C Sort out the chosen points C if ((.not. quit1) .and. (.not. quit2)) then xmin = min(xid(1),xid(2)) xmax = max(xid(1),xid(2)) ymin = min(yid(1),yid(2)) ymax = max(yid(1),yid(2)) else status = 1 end if end subroutine reject(ncomax,nco,nx,ny,picture,xaxis,yaxis,params, : errpar,mon,resmin,resmax) C+-------------------------------------------------------------------------------------- C REJECT --- Reject bulk fits C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) NCOMAX: (integer) Maximum number of coefs C (>) NCO: (integer) Number of coefs C (>) NX: (integer) X dimension of the data C (>) NY: (integer) Y dimension of the data C (>) PICTURE: (logical) TRUE if a 2-d plot is on the screen C (>) XAXIS: (real array XAXIS(NX)) X axis data C (>) YAXIS: (real array YAXIS(NY)) Y axis data C (!) PARAMS: (real array PARAMS(NCOMAX+4,NX,NY)) Final output array C (!) ERRPAR: (real array ERRPAR(NCOMAX+4,NX,NY)) Output error array C (!) MON: (integer array MON(4,NX,NY)) Monitor array C (>) RESMIN: (real array RESMIN(NCOMAX)) Minimum values for params C (>) RESMAX: (real array RESMAX(NCOMAX)) Maximum values for params C C C AUTHOR: Jim Lewis -- RGO C DATE: 28-NOV-1990 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer ncomax,nco,nx,ny integer mon(4,nx,ny) real params(ncomax+4,nx,ny),errpar(ncomax+4,nx,ny),xaxis(nx), : yaxis(ny),resmin(ncomax),resmax(ncomax) logical picture C C Functions C integer ich_len logical par_batch logical par_quest C C Local variables C integer status,i,j,k,index,nreject real xmin,xmax,ymin,ymax,x,y,vmin(2),vmax(2),start(2),end(2), : dummy(2),magic_float logical use_curs,inside,in,stop,rej_dom(2),ditch,batch character message*80 C C Parameters for DSA and PAR routines... C character range(4)*7,units(4)*7 data range /'rng_amp','rng_pos','rng_sig','rng_bac'/ data units /'counts ','microns','microns','counts'/ C C Numeric ranges C include 'NDP_NUMERIC_RANGES' C C Initialise the domain array C do i = 1,2 rej_dom(i) = .false. dummy(i) = 0.0 end do call dsa_get_flag_value('FLOAT',magic_float,status) if (status .ne. 0) go to 500 C C Now prompt for the domains to be considered C call dsa_wruser(' ') call dsa_wrflush call dsa_wruser : ('You may reject the fits in the following domains:') call dsa_wrflush call dsa_wruser('(1) The x-y plane') call dsa_wrflush call dsa_wruser('(2) By the results') call dsa_wrflush batch = par_batch() if (.not. batch) call par_cnpar('rej_dom') call par_rdary('rej_dom',0.0,1.0,'n',' ',2,2,dummy) do i = 1,2 index = nint(dummy(i)) if (index .eq. 1) rej_dom(i) = .true. end do nreject = 0 C C If you want to reject in xy space, then do so... C if (rej_dom(1)) then C C Define regions with cursor? C if (picture) then use_curs = par_quest('Define region with cursor?',.true.) else use_curs = .false. end if C C Do you zap out the bits inside or outside the frame? C inside = par_quest('Delete the bit INSIDE the square?',.true.) C C Define regions to reject until you want to quit C stop = .false. status = 0 do while (.not. stop) C C If you're using a 2-d plot, try and get the limits from the cursor C if (use_curs) then call dsa_wruser(' ') call dsa_wrflush call dsa_wruser('Define two points with the cursor, ') call dsa_wruser('"Q" to quit') call dsa_wrflush call curse_limits(nx,ny,xaxis,yaxis,xmin,xmax, : ymin,ymax,status) if (status .ne. 0) stop = .true. if ((xmin .lt. 1.0) .or. (xmax .gt. float(nx)) .or. : (ymin .lt. 1.0) .or. (ymax .gt. float(ny))) stop = .true. end if C C If you're not using a 2-d plot, or if you are and not using the C cursor to define your region, then input the limits by hand. C if (((.not. picture) .or. ((picture) : .and. (.not. use_curs))) .and. (.not. stop)) then vmin(1) = float(min_int) vmin(2) = float(min_int) vmax(1) = float(max_int) vmax(2) = float(max_int) if (.not. batch) then call par_cnpar('start') call par_cnpar('end') end if call ndp_par_rdary('start',vmin,vmax,' ','pixels',2,2, : start) call ndp_par_rdary('end',vmin,vmax,' ','pixels',2,2,end) xmin = start(1) ymin = start(2) xmax = end(1) ymax = end(2) if ((xmin .lt. 1.0) .or. (xmax .gt. float(nx)) .or. : (ymin .lt. 1.0) .or. (ymax .gt. float(ny))) stop = .true. end if C C Draw a box of the defined limits C if (.not. stop) then if (picture) then call pgsci(2) call pgscr(2,1.0,0.0,0.0) call pgwindow(1.0,float(nx),1.0,float(ny)) call pgmove(xmin,ymin) call pgdraw(xmax,ymin) call pgdraw(xmax,ymax) call pgdraw(xmin,ymax) call pgdraw(xmin,ymin) call pgsci(1) end if C C Now zap out the bits you don't want... C do j = 1,ny y = float(j) do i = 1,nx x = float(i) in = .false. if ((x .gt. xmin) .and. (x .lt. xmax) .and. : (y .gt. ymin) .and. (y .lt. ymax)) in = .true. if (((in) .and. (inside)) .or. ((.not. in) .and. : (.not. inside))) then do k = 1,ncomax params(k,i,j) = magic_float errpar(k,i,j) = magic_float end do if (mon(3,i,j) .gt. 0) then mon(3,i,j) = -1 nreject = nreject + 1 end if end if end do end do end if end do end if C C If you wish to reject by the results, then do so... C if (rej_dom(2)) then C C Read the ranges to be considered now... C do i = 1,nco if (.not. batch) call par_cnpar(range(i)) call par_rdary(range(i),min_float,max_float,'a',units(i),2, : 2,dummy) resmin(i) = dummy(1) resmax(i) = dummy(2) end do C C Now go through and zap the bad ones... C do j = 1,ny do i = 1,nx if ((mon(3,i,j) .gt. 0) .and. : ((mon(1,i,j) .gt. 0) .or. (mon(1,i,j) .eq. -1))) then ditch = .false. do k = 1,nco if ((params(k,i,j) .lt. resmin(k)) .or. : (params(k,i,j) .gt. resmax(k))) ditch = .true. end do if (ditch) then do k = 1,nco params(k,i,j) = magic_float errpar(k,i,j) = magic_float end do mon(3,i,j) = -1 nreject = nreject + 1 end if end if end do end do end if C C Now tell how many were zapped this time for one reason or another... C write(message,100) nreject 100 format(i8,' spectra were rejected') call dsa_wruser(message(:ich_len(message))) call dsa_wrflush 500 end subroutine edit(nx,ny,nz,nco,ncomax,picture,zaxis,params,errpar, : data,err,niter,nkeep,fit_type,how_guess,alam0,dlam, : mon,used,xmin,xmax,ymin,ymax,work1,work2,iwt,boxsize) C+-------------------------------------------------------------------------------------- C EDIT --- Edit the fits by hand C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) NX: (integer) X dimension of the data C (>) NY: (integer) Y dimension of the data C (>) NZ: (integer) Z dimension of the data C (>) NCO: (integer) Number of coefs C (>) NCOMAX: (integer) Maximum number of coefs C (>) PICTURE: (logical) TRUE if a 2-d plot is on the screen C (>) ZAXIS: (real array ZAXIS(NZ)) Spectral (X) axis values C (!) PARAMS: (real array PARAMS(NCOMAX+4,NX,NY)) Final output array C (!) ERRPAR: (real array ERRPAR(NCOMAX+4,NX,NY)) Output error array C (>) DATA: (real array DATA(NZ,NKEEP)) Input data array C (>) ERR: (real array ERR(NZ,NKEEP)) Input error array C (>) NITER: (integer) Maximum number of iterations C (>) NKEEP: (integer) Number of spectra kept C (>) FIT_TYPE (integer) Flag for fitting function C (>) HOW_GUESS (integer) Method for guessing C (>) ALAM0: (real) Wavelength zeropoint C (>) DLAM: (real) Wavelength increment C (!) MON: (integer array MON(4,NX,NY)) Monitor array C (>) USED: (integer array USED(3,NKEEP)) Flags for used spectra C (>) XMIN: (real) Lower limit on X to be considered C (>) XMAX: (real) Upper limit on X to be considered C (>) YMIN: (real) Lower limit on Y to be considered C (>) YMAX: (real) Upper limit on Y to be considered C (W) WORK1: (real array WORK1(NZ)) Workspace for plotting routines C (W) WORK2: (real array WORK2(NZ)) Workspace for plotting routines C (>) IWT: (integer) Flag for weighting scheme C (>) BOXSIZE: (integer) Size of binning box C C C AUTHOR: Jim Lewis -- RGO C DATE: 28-NOV-1990 C UPDATE: 19-DEC-1990 -- Added options I,C,S. (JRL) C 18-JAN-1991 -- Added option W (JRL) C 08-SEP-1992 -- Added option L (JRL) C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nx,ny,nz,nkeep,ncomax,nco,niter,how_guess,fit_type, : iwt,boxsize integer mon(4,nx,ny),used(3,nkeep) real params(ncomax+4,nx,ny),errpar(ncomax+4,nx,ny), : zaxis(nz),work1(nz),work2(nz),alam0,dlam,data(nz,nkeep), : err(nz,nkeep),xmin,xmax,ymin,ymax logical picture C C Functions C integer dsa_typesize integer dyn_element integer ich_fold integer ich_len logical par_quest C C Local variables C integer np,status,totspace,wpt1,icount,jx,jy,guess_status, : address,slot,dumint,nplotted,bytes_int,bytes_float,ix,iy, : nspec,npage,slot2,cptr,eptr,zptr,coptr,alptr,lptr,nfail,nrun, : icx,icy,unit,nlast,ntriv_min,wpt2,slot3,length real deltax,deltay,margin,xvstart,xdev,ydev,yvstart,zmin,zmax,xx, : yy,xv1,xv2,yv1,yv2,dmin,dmax,fsr_wave,xv1new,xv2new,yv1new, : yv2new,fsr,dchi_max logical quit,stop character option*1,message*80,helpfile*80,harddev*16,pgstate*20 C C Common block with FSR C common /free/ fsr C C Common block with convergence parameters C common /conv/ dchi_max,ntriv_min C C Common block for dynamic memory. C include 'DYNAMIC_MEMORY' C C Initialise workspace slots and some workspace variables C slot = 0 slot2 = 0 bytes_float = dsa_typesize('float',status) bytes_int = dsa_typesize('int',status) C C Set up screen C call clear_screen(status) if (status .ne. 0) go to 500 picture = .false. call pre_plot(nz,zaxis,zmin,zmax,alam0,dlam,nkeep,xvstart,yvstart, : margin,deltax,deltay,xdev,ydev,np,slot,wpt1,status) if (status .ne. 0) go to 500 call plot_panel(np,xvstart,yvstart,deltax,deltay,margin,ncomax, : params,mon,data,nkeep,nx,ny,nz,xmin,xmax,ymin,ymax,zmin,zmax, : zaxis,nplotted,nlast,dlam,work1,work2,dynamic_mem(wpt1),xv1,xv2, : yv1,yv2,.true.,.false.,.false.) C C FSR in wavelength units C fsr_wave = fsr*dlam C C Count the spectra C icount = 0 C C Keep at it until you've edited them all or until you're too tired C to continue... C stop = .false. do while ((icount .lt. nkeep) .and. (.not. stop)) quit = .false. C C First, plot up some spectra and their fits... C call plot_panel(np,xvstart,yvstart,deltax,deltay,margin,ncomax, : params,mon,data,nkeep,nx,ny,nz,xmin,xmax,ymin,ymax,zmin,zmax, : zaxis,nplotted,nlast,dlam,work1,work2,dynamic_mem(wpt1),xv1, : xv2,yv1,yv2,.false.,.false.,.false.) icount = icount + nplotted C C Now select the spectra you want to deal with individually C call dsa_wruser('Select individual spectra with cursor') call dsa_wrflush do while (.not. quit) call pgvsize(xv1,xv2,yv1,yv2) call minmax(nz,data(1,nlast),dmin,dmax) call pgwindow(zmin,zmax,dmin,dmax) call pgcurse(xx,yy,option) dumint = ich_fold(option) C C Option to move to the next page (Next) C if (option .eq. 'N') then quit = .true. C C Option to delete a fit (Zap) C else if (option .eq. 'Z') then C C First find out which spectrum it is... C call find_cell(np,nz,xx,yy,dynamic_mem(wpt1), : xv1,yv1,deltax,deltay,margin,zmin,zmax,data(1,nlast), : ix,iy,npage,nspec) C C Set up the plot device so you can erase the fit... C if (nspec .gt. 0) then icx = used(1,nspec) icy = used(2,nspec) if (mon(3,icx,icy) .gt. 0) then xv1new = margin + float(ix-1)*(deltax + 2.0*margin) yv1new = margin + float(iy-1)*(deltay + 2.0*margin) xv2new = xv1new + deltax yv2new = yv1new + deltay call which_spectrum(dynamic_mem(wpt1),np,npage,jx, : jy) call pgvsize(xv1new,xv2new,yv1new,yv2new) call minmax(nz,data(1,nspec),dmin,dmax) call pgwindow(zmin,zmax,dmin,dmax) if (fit_type .eq. 1) then call gauss_eval(params(1,jx,jy),ncomax,zaxis,nz, : fsr_wave,work1) else if (fit_type .eq. 2) then call cauchy_eval(params(1,jx,jy),ncomax,zaxis, : nz,fsr_wave,work1) end if C C Now erase the fit C call pgsci(0) call pgline(nz,zaxis,work1) call pgsci(1) call zap_fit(nx,ny,np,npage,ncomax, : dynamic_mem(wpt1),fit_type,mon,params,errpar) end if else call dsa_wruser('There is no spectrum here!!') call dsa_wrflush end if C C Option to fit the spectrum (Fit) C else if (option .eq. 'F') then C C First get some workspace C totspace = (2*nco + nz + : 2*nco*nco)*bytes_float + nco*bytes_int call dsa_get_workspace(totspace,address,slot2,status) if (status .ne. 0) then call dsa_wruser : ('Couldn''t allocate workspace in editing routine') call dsa_wrflush go to 500 end if cptr = dyn_element(address) eptr = cptr + nco*bytes_float zptr = eptr + nco*bytes_float coptr = zptr + nz*bytes_float alptr = coptr + nco*nco*bytes_float lptr = alptr + nco*nco*bytes_float C C Find out which cell on the page you're dealing with C call find_cell(np,nz,xx,yy,dynamic_mem(wpt1), : xv1,yv1,deltax,deltay,margin,zmin,zmax,data(1,nlast), : ix,iy,npage,nspec) if (nspec .gt. 0) then C C If you are guessing by choosing another spectrum, then do that C here C if (how_guess .eq. 5) then guess_status = 0 call guess_other(nx,ny,nz,np,dynamic_mem(wpt1),xv1, : yv1,deltax,deltay,margin,zmin,zmax,data(1,nlast), : nco,ncomax,nkeep,mon,used,params,alam0,dlam, : dynamic_mem(cptr),guess_status) if (guess_status .ne. 0) go to 450 end if C C Reset the viewport for the spectrum of interest C xv1 = margin + float(ix-1)*(deltax + 2.0*margin) yv1 = margin + float(iy-1)*(deltay + 2.0*margin) xv2 = xv1 + deltax yv2 = yv1 + deltay call pgvsize(xv1,xv2,yv1,yv2) C C Fit the spectrum C nfail = 0 nrun = 0 call fit_spec(nco,ncomax,dynamic_mem(lptr),nfail,nrun, : nx,ny,nz,how_guess,nkeep,nspec,data(1,nspec), : err(1,nspec),used,dynamic_mem(zptr),fsr, : dynamic_mem(cptr),dynamic_mem(eptr),niter, : dynamic_mem(coptr),dynamic_mem(alptr),fit_type,alam0, : dlam,params) C C Store the fit if it succeded and issue appropriate warning messages C Draw the fit, if it has been accepted C if (nfail .ne. 0) then write(message,100) niter 100 format('Spectrum failed to converge within',i4, : ' iterations') call dsa_wruser(message(:ich_len(message))) call dsa_wrflush call zap_fit(nx,ny,np,npage,ncomax, : dynamic_mem(wpt1),fit_type,mon,params,errpar) else if (nrun .ne. 0) then call dsa_wruser : ('Spectrum had a runaway solution') call dsa_wrflush end if call accept_fit(nx,ny,np,npage,dynamic_mem(wpt1), : dynamic_mem(cptr),dynamic_mem(eptr),ncomax,nco, : fit_type,boxsize,mon,params,errpar) call minmax(nz,data(1,nspec),dmin,dmax) call pgwindow(zmin,zmax,dmin,dmax) if (fit_type .eq. 1) then call gauss_eval(dynamic_mem(cptr),nco,zaxis,nz, : fsr_wave,work1) else if (fit_type .eq. 2) then call cauchy_eval(dynamic_mem(cptr),nco,zaxis,nz, : fsr_wave,work1) end if C C Now plot the fit as a green line... C call pgsci(3) call pgscr(3,0.0,1.0,0.0) call pgline(nz,zaxis,work1) call pgsci(1) end if call dsa_free_workspace(slot2,status) slot2 = 0 else call dsa_wruser('There is no spectrum here!!') call dsa_wrflush end if C C Option to change the fitting function (Type) C else if (option .eq. 'T') then call type_fit(.true.,fit_type,nco) C C Option to change the guessing algorithm (Guess) C else if (option .eq. 'G') then call find_guess(how_guess,.false.,.true.) C C Option to redraw the current panel in the event you have schmucked C it up...(Redraw) C else if (option .eq. 'R') then call plot_panel(np,xvstart,yvstart,deltax,deltay,margin, : ncomax,params,mon,data,nkeep,nx,ny,nz,xmin,xmax,ymin, : ymax,zmin,zmax,zaxis,nplotted,nlast,dlam,work1,work2, : dynamic_mem(wpt1),xv1,xv2,yv1,yv2,.false.,.true.,.false.) C C Option to move to the previous panel (Previous) C else if (option .eq. 'P') then icount = icount - nplotted call plot_panel(np,xvstart,yvstart,deltax,deltay,margin, : ncomax,params,mon,data,nkeep,nx,ny,nz,xmin,xmax,ymin, : ymax,zmin,zmax,zaxis,nplotted,nlast,dlam,work1,work2, : dynamic_mem(wpt1),xv1,xv2,yv1,yv2,.false.,.false.,.true.) C C Print up a help menu (Help) C else if (option .eq. 'H') then call dsa_open_text_file('TAU_FITS_EDIT','.TXT','old', : .false.,unit,helpfile,status) if (status .ne. 0) then status = 0 call dsa_wruser : ('Sorry, I can''t find the HELP file') call dsa_wrflush else call dsa_wruser : ('****The options available in the menu are:') call dsa_wrflush call dsa_wruser(' ') call dsa_wrflush read(unit,105,iostat=status) message 105 format(a80) do while (status .eq. 0) call dsa_wruser(message(:ich_len(message))) call dsa_wrflush read(unit,105,iostat=status) message end do call dsa_wruser(' ') call dsa_wrflush status = 0 close(unit) call dsa_free_lu(unit,status) end if C C Option to print out setup info (Info) C else if (option .eq. 'I') then call info(fit_type,iwt,how_guess,boxsize,niter, : ntriv_min,dchi_max) C C Option to change convergence parameters (Convergence) C else if (option .eq. 'C') then call conv_pars(ntriv_min,dchi_max,niter) C C Option to show the fit parameters for a given fit (Show) C else if (option .eq. 'S') then call find_cell(np,nz,xx,yy,dynamic_mem(wpt1), : xv1,yv1,deltax,deltay,margin,zmin,zmax,data(1,nlast), : ix,iy,npage,nspec) if (nspec .gt. 0) then icx = used(1,nspec) icy = used(2,nspec) write(message,110) icx,icy 110 format('Pixel (',i3,',',i3,')') call dsa_wruser(message(:ich_len(message))) call dsa_wrflush totspace = 4*ncomax*bytes_float call dsa_get_workspace(totspace,address,slot3,status) if (status .ne. 0) then status = 0 go to 450 end if wpt2 = dyn_element(address) call gen_fill(totspace,0,dynamic_mem(wpt2)) call show(nx,ny,float(icx),float(icx),float(icy), : float(icy),params,errpar,ncomax,mon,dynamic_mem(wpt2)) call dsa_free_workspace(slot3,status) if (status .ne. 0) then status = 0 go to 450 end if else call dsa_wruser('There is no spectrum here!!') call dsa_wrflush end if C C Option to delete all fits from a panel (Wipe) C else if (option .eq. 'W') then call wipe_out(nx,ny,np,ncomax,dynamic_mem(wpt1),fit_type, : mon,params,errpar) call plot_panel(np,xvstart,yvstart,deltax,deltay,margin, : ncomax,params,mon,data,nkeep,nx,ny,nz,xmin,xmax,ymin, : ymax,zmin,zmax,zaxis,nplotted,nlast,dlam,work1,work2, : dynamic_mem(wpt1),xv1,xv2,yv1,yv2,.false.,.true.,.false.) C C Option to get a hard copy (Laser) C else if (option .eq. 'L') then C C Close current plot device and open hardcopy device C call pgqinf('state',pgstate,length) if (pgstate .eq. 'OPEN') call pgend call var_getchr('harddev',0,0,harddev,status) if (status .ne. 0) then call dsa_wruser : ('No hardcopy device selected -- use HARD command') call dsa_wrflush go to 500 end if call pgbegin(0,harddev,1,1) C C Do the plot C call plot_panel(np,xvstart,yvstart,deltax,deltay,margin, : ncomax,params,mon,data,nkeep,nx,ny,nz,xmin,xmax,ymin, : ymax,zmin,zmax,zaxis,nplotted,nlast,dlam,work1,work2, : dynamic_mem(wpt1),xv1,xv2,yv1,yv2,.false.,.true.,.false.) C C Message to the user C call dsa_wruser('Hardcopy done. Now please wait a sec ') call dsa_wruser(' while I reopen the soft device') call dsa_wrflush C C Clear the screen and replot the data in case the user wants C to fiddle with this panel again C call clear_screen(status) if (status .ne. 0) go to 500 call plot_panel(np,xvstart,yvstart,deltax,deltay,margin, : ncomax,params,mon,data,nkeep,nx,ny,nz,xmin,xmax,ymin, : ymax,zmin,zmax,zaxis,nplotted,nlast,dlam,work1,work2, : dynamic_mem(wpt1),xv1,xv2,yv1,yv2,.false.,.true.,.false.) C C Option to get out of EDIT altogether (Quit) C else if (option .eq. 'Q') then quit = par_quest('Do you REALLY want to quit?',.true.) if (quit) stop = .true. end if 450 end do end do C C Tidy and exit C 500 if (slot .ne. 0) call dsa_free_workspace(slot,status) if (slot2 .ne. 0) call dsa_free_workspace(slot2,status) end subroutine pre_plot(nz,zaxis,zmin,zmax,alam0,dlam,nkeep,xvstart, : yvstart,margin,deltax,deltay,xdev,ydev,np,slot,wpt1,status) C+-------------------------------------------------------------------------------------- C PRE_PLOT --- Set up the screen for panel of spectral plots C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) NZ: (integer) Z dimension of the data C (>) ZAXIS: (real array ZAXIS(NZ)) Z axis data array C (>) ZMIN: (real) Minimum value of Z C (>) ZMAX: (real) Maximum value of Z C (>) ALAM0: (real) Wavelength zeropoint C (>) DLAM: (real) Wavelength increment C (>) NKEEP: (integer) Total number of spectra analysed C (<) XVSTART: (real) X pos on viewport for next plot C (<) YVSTART: (real) Y pos on viewport for next plot C (<) MARGIN: (real) Size of margin for each plot C (<) DELTAX: (real) X size of plot C (<) DELTAY: (real) Y size of plot C (<) XDEV: (real) X size of view surface C (<) YDEV: (real) Y size of view surface C (<) NP: (integer) Number of plots on a side C (<) SLOT: (integer) Workspace slot C (<) WPT1: (integer) Workspace pointer C (<) STATUS: (integer) Status variable C C C AUTHOR: Jim Lewis -- RGO C DATE: 28-NOV-1990 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer status,np,slot,wpt1,nz,nkeep real deltax,deltay,xvstart,yvstart,margin,zaxis(nz),zmin,zmax, : alam0,dlam,xdev,ydev C C Functions C integer dsa_typesize integer dyn_element C C Local variables C integer ignore,ntot,bytes_int,totspace,address,i real dummy2,dummy,x1,x2,y1,y2 C C Common block for dynamic memory. C include 'DYNAMIC_MEMORY' C C Get number of plots on a side C call var_getnum('num_plots',0,0,dummy2,ignore) call par_qnum('Number of plots on a side',1.0,10.0,dummy2,.true., : ' ',dummy) call var_setnum('num_plots',0,0,dummy,ignore) np = nint(dummy) ntot = np**2 C C Get some information about the plot device and set up C plot sizes and margins C call pgvstand call pgqvp(1,x1,x2,y1,y2) xdev = x2 + x1 ydev = y2 + y1 deltax = xdev/(float(np)*1.2) margin = 0.1*deltax deltay = ydev/float(np) - 2.0*margin xvstart = margin yvstart = margin C C Get some workspace and initialise it C bytes_int = dsa_typesize('int',status) totspace = 3*ntot*bytes_int call dsa_get_workspace(totspace,address,slot,status) if (status .ne. 0) then call dsa_wruser : ('Couldn''t allocate workspace in pre-plotting routine') call dsa_wrflush go to 500 end if wpt1 = dyn_element(address) call gen_fill(totspace,0,dynamic_mem(wpt1)) C C Set up spectral array C do i = 1,nz zaxis(i) = alam0 + float(i)*dlam end do zmin = min(zaxis(1),zaxis(nz)) zmax = max(zaxis(1),zaxis(nz)) 500 end subroutine which_spectrum(monitor,np,npage,jx,jy) C+-------------------------------------------------------------------------------------- C WHICH_SPECTRUM --- Find which spectrum has been indicated with cursor C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) MONITOR: (integer array MONITOR(3,NP*NP)) Plot monitor array C (>) NP: (integer) Number of plots on a side C (>) NPAGE: (integer) Number of the indicated spectrum on the page C (<) JX: (integer) X position of spectrum C (<) JY: (integer) Y position of spectrum C C C AUTHOR: Jim Lewis -- RGO C DATE: 28-NOV-1990 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer np,npage,jx,jy integer monitor(3,np*np) C C Get the spectrum position C jx = monitor(1,npage) jy = monitor(2,npage) end subroutine accept_fit(nx,ny,np,npage,monitor,coefs,errcoefs, : ncomax,nco,fit_type,boxsize,mon,params,errpar) C+-------------------------------------------------------------------------------------- C ACCEPT_FIT --- Accept an individual fit done by hand C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) NX: (integer) X dimension of the data C (>) NY: (integer) Y dimension of the data C (>) NP: (integer) Number of plots on a side C (>) NPAGE: (integer) Number of the indicated spectrum on the page C (>) MONITOR: (integer array MONITOR(3,NP*NP)) Plot monitor array C (>) COEFS: (real array COEFS(NCO)) Fit coefs C (>) ERRCOEFS: (real array ERRCOEFS(NCO)) Error in fit coefs C (>) NCOMAX: (integer) Maximum number of coefs C (>) NCO: (integer) Number of coefs C (>) FIT_TYPE: (integer) Flag for fitting function C (>) BOXSIZE: (integer) Size of binning box C (!) MON: (integer array MON(4,NX,NY)) Monitor array C (<) PARAMS: (real array PARAMS(NCOMAX+4,NX,NY)) Final output array C (<) ERRPAR: (real array ERRPAR(NCOMAX+4,NX,NY)) Output array errors C C C AUTHOR: Jim Lewis -- RGO C DATE: 28-NOV-1990 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nx,ny,npage,np,ncomax,nco,fit_type,boxsize integer monitor(3,np*np),mon(4,nx,ny) real coefs(nco),errcoefs(nco),params(ncomax+4,nx,ny), : errpar(ncomax+4,nx,ny) C C Local variables C integer i,ix,iy C C Which spectrum? C call which_spectrum(monitor,np,npage,ix,iy) C C Write the fit coefs to the output array C do i = 1,nco params(i,ix,iy) = coefs(i) errpar(i,ix,iy) = errcoefs(i) end do mon(2,ix,iy) = fit_type mon(3,ix,iy) = 2 mon(4,ix,iy) = boxsize end subroutine zap_fit(nx,ny,np,npage,ncomax,monitor,fit_type, : mon,params,errpar) C+-------------------------------------------------------------------------------------- C ZAP_FIT --- Zap an individual fit done by hand C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) NX: (integer) X dimension of the data C (>) NY: (integer) Y dimension of the data C (>) NP: (integer) Number of plots on a side C (>) NPAGE: (integer) Number of the indicated spectrum on the page C (>) NCOMAX: (integer) Maximum number of coefs C (>) MONITOR: (integer array MONITOR(3,NP*NP)) Plot monitor array C (>) FIT_TYPE: (integer) Flag for fitting function C (>) MON: (integer array MON(4,NX,NY)) Monitor array C (<) PARAMS: (real array PARAMS(NCOMAX+4,NX,NY)) Final output array C (<) ERRPAR: (real array ERRPAR(NCOMAX+4,NX,NY)) Output array errors C C C AUTHOR: Jim Lewis -- RGO C DATE: 28-NOV-1990 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nx,ny,npage,ncomax,np,fit_type integer monitor(3,np*np),mon(4,nx,ny) real params(ncomax+4,nx,ny),errpar(ncomax+4,nx,ny) C C Local variables C integer i,ix,iy,status real magic_float C C Magic value C call dsa_get_flag_value('FLOAT',magic_float,status) if (status .ne. 0) go to 500 C C Which spectrum? C call which_spectrum(monitor,np,npage,ix,iy) C C Replace the deleted spectrum parameters with magic values C do i = 1,ncomax params(i,ix,iy) = magic_float errpar(i,ix,iy) = magic_float end do mon(2,ix,iy) = fit_type mon(3,ix,iy) = -2 500 end subroutine find_cell(np,nz,xx,yy,monitor,xv1,yv1, : deltax,deltay,margin,zmin,zmax,data,ix,iy,npage,nspec) C+-------------------------------------------------------------------------------------- C FIND_CELL --- Find the cell on the page at which the cursor is pointing C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) NP: (integer) Number of plots on a side C (>) NZ: (integer) Z dimension of the data C (>) XX: (real) Cursor X position C (>) YY: (real) Cursor Y position C (>) MONITOR: (integer array MONITOR(3,NP*NP)) Plot monitor array C (>) XV1: (real) X position of lower left hand corner of screen C (>) YV1: (real) Y position of lower left hand corner of screen C (>) DELTAX: (real) X size of plot in inches C (>) DELTAY: (real) Y size of plot in inches C (>) MARGIN: (integer) Size of margin in inches C (>) ZMIN: (real) Minimum value of Z C (>) ZMAX: (real) Maximum value of Z C (>) DATA: (real array DATA(NZ)) Data array of last spectrum plotted C (<) IX: (integer) X cell position on page of plot C (<) IY: (integer) Y cell position on page of plot C (<) NPAGE: (integer) Number of the cell on the page C (<) NSPEC: (integer) Number of the spectrum C C C AUTHOR: Jim Lewis -- RGO C DATE: 28-NOV-1990 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer np,ix,iy,nspec,npage,nz integer monitor(3,np*np) real xx,yy,xv1,yv1,deltax,deltay,zmin,zmax,data(nz),margin C C Local variables C real dmin,dmax,xscale,yscale,xp,yp C C First find the scale in inches/world coords of the last graph plotted C call minmax(nz,data,dmin,dmax) xscale = deltax/(zmax-zmin) yscale = deltay/(dmax-dmin) C C Now find the position of the cursor in inches from the lower left hand C side of the view surface C xp = xv1 + (xx - zmin)*xscale yp = yv1 + (yy - dmin)*yscale C C Now find the cell position of the cursor C ix = int(xp/(deltax + 2.0*margin)) + 1 iy = int(yp/(deltay + 2.0*margin)) + 1 npage = np*(iy-1) + ix nspec = monitor(3,npage) end subroutine wipe_out(nx,ny,np,ncomax,monitor,fit_type, : mon,params,errpar) C+-------------------------------------------------------------------------------------- C WIPE_OUT --- Zap a whole panel of fits C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) NX: (integer) X dimension of the data C (>) NY: (integer) Y dimension of the data C (>) NP: (integer) Number of plots on a side C (>) NCOMAX: (integer) Maximum number of coefs C (>) MONITOR: (integer array MONITOR(3,NP*NP)) Plot monitor array C (>) FIT_TYPE: (integer) Flag for fitting function C (>) MON: (integer array MON(4,NX,NY)) Monitor array C (<) PARAMS: (real array PARAMS(NCOMAX+4,NX,NY)) Final output array C (<) ERRPAR: (real array ERRPAR(NCOMAX+4,NX,NY)) Output array errors C C C AUTHOR: Jim Lewis -- RGO C DATE: 18-JAN-1991 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nx,ny,np,ncomax,fit_type integer monitor(3,np*np),mon(4,nx,ny) real params(ncomax+4,nx,ny),errpar(ncomax+4,nx,ny) C C Local variables C integer npage,ix,iy C C Loop through and zap any spectra which have a valid fit... C do npage = 1,np*np call which_spectrum(monitor,np,npage,ix,iy) if ((ix .ne. 0) .and. (iy .ne. 0)) then if (mon(3,ix,iy) .gt. 0) then call zap_fit(nx,ny,np,npage,ncomax,monitor,fit_type, : mon,params,errpar) end if end if end do end subroutine fit_spec(nco,ncomax,lista,nfail,nrun,nx,ny,nz, : how_guess,nkeep,nspec,data,err,used,lamda,fsr,coefs,errcoefs, : niter,covar,alpha,fit_type,alam0,dlam,params) C+-------------------------------------------------------------------------------------- C FIT_SPEC --- Routine to fit an individual spectrum C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) NCO: (integer) Number of coefs C (>) NCOMAX: (integer) Maximum number of coefs C (>) LISTA: (integer array LISTA(NCO)) List of coefs to be varied. C (!) NFAIL: (integer) Number of fits which have failed C (!) NRUN: (integer) Number of fits which have runaway solutions C (>) NX: (integer) X dimension of the data C (>) NY: (integer) Y dimension of the data C (>) NZ: (integer) Z dimension of the data C (>) HOW_GUESS: (integer) Guessing algorithm C (>) NKEEP: (integer) Number of spectra kept C (>) NSPEC: (integer) Number of the spectrum C (>) DATA: (real array DATA(NZ)) Input data array C (>) ERR: (real array ERR(NZ)) Input error array C (>) USED: (integer array USED(3,NKEEP)) Flags for used spectra C (>) LAMDA: (real array LAMDA(NZ)) Spectral axis data C (>) FSR: (real) Free spectral range C (<) COEFS: (real array COEFS(NCO)) Fitted coefs C (<) ERRCOEFS: (real array ERRCOEFS(NCO)) Error in fitted coefs C (>) NITER: (integer) Maximum number of iterations for fit. C (W) COVAR: (real array COVAR(NCO,NCO)) Covariance matrix workspace C (W) ALPHA: (real array ALPHA(NCO,NCO)) Curvature matrix workspace C (>) FIT_TYPE: (integer) Flag for fitting function C (>) ALAM0: (real) Wavelength zeropoint C (>) DLAM: (real) Wavelength increment C (>) PARAMS: (real array PARAMS(NCOMAX+4,NX,NY)) Final output array C C C AUTHOR: Jim Lewis -- RGO C DATE: 28-NOV-1990 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nco,nfail,nrun,nx,ny,nz,how_guess,nkeep,nspec,niter, : ncomax,fit_type integer lista(nco),used(3,nkeep) real data(nz),err(nz),lamda(nz),fsr,coefs(nco), : covar(nco,nco),alpha(nco,nco),alam0,dlam,errcoefs(nco), : params(ncomax+4,nx,ny) C C Functions C integer dsa_typesize integer dyn_element C C Local variables C integer i,bytes_float,bytes_int,totspace,status,slot,iwpt1, : iwpt2,iwpt3,iwpt4,iwpt5,iwpt6,iwpt7,address,slot2,wptr logical quick C C Common block for dynamic memory. C include 'DYNAMIC_MEMORY' C C A few initialisations for the fitting routine... C do i = 1,nco lista(i) = i end do status = 0 slot = 0 slot2 = 0 C C Get workspace for fitting routines... C bytes_float = dsa_typesize('float',status) bytes_int = dsa_typesize('int',status) totspace = 4*nco*bytes_float + 3*nco*bytes_int : + 2*nz*bytes_float call dsa_get_workspace(totspace,address,slot,status) if (status .ne. 0) go to 500 iwpt1 = dyn_element(address) iwpt2 = iwpt1 + nco*bytes_float iwpt3 = iwpt2 + nco*bytes_float iwpt4 = iwpt3 + nco*bytes_float iwpt5 = iwpt4 + nco*bytes_float iwpt6 = iwpt5 + 3*nco*bytes_int iwpt7 = iwpt6 + nz*bytes_float C C Now do the fit -- first take a guess at the coefs C if (how_guess .eq. 1) then call guess_moms(data,nz,fsr,coefs) else if (how_guess .eq. 2) then call guess_previous(data,nx,ny,nz,fsr,nco,ncomax, : nspec,nkeep,used,params,alam0,dlam,coefs) else if ((how_guess .eq. 3) .or. (how_guess .eq. 4)) then totspace = nz*bytes_float call dsa_get_workspace(totspace,address,slot2,status) if (status .ne. 0) go to 500 wptr = dyn_element(address) quick = .true. if (how_guess .eq. 3) quick = .false. call guess_by_hand(quick,data,nz,lamda,dynamic_mem(wptr),coefs, : status) call dsa_free_workspace(slot2,status) slot2 = 0 end if C C Check to see if the line is near an edge C call fiddle(coefs,data,err,nz,fsr,dynamic_mem(iwpt6), : dynamic_mem(iwpt7),lamda,status) C C Now do the fit C call sub_fit(niter,nz,lamda,dynamic_mem(iwpt6), : dynamic_mem(iwpt7),nco,lista,covar,alpha,fit_type,alam0, : dlam,dynamic_mem(iwpt1),dynamic_mem(iwpt2), : dynamic_mem(iwpt3),dynamic_mem(iwpt4),dynamic_mem(iwpt5), : coefs,errcoefs,nfail,nrun,status) C C Tidy and exit C 500 status = 0 if (slot .ne. 0) call dsa_free_workspace(slot,status) if (slot2 .ne. 0) call dsa_free_workspace(slot2,status) end subroutine plot_panel(np,xvstart,yvstart,deltax,deltay,margin, : ncomax,params,mon,data,nkeep,nx,ny,nz,xmin,xmax,ymin,ymax,zmin, : zmax,zaxis,nplotted,nlast,dlam,work1,work2,monitor,xv1,xv2, : yv1,yv2,setup,redraw,backwards) C+-------------------------------------------------------------------------------------- C PLOT_PANEL --- Plot a 2-d panel of spectra C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) NP: (integer) Number of plots on a side C (>) XVSTART: (real) Where in the viewport to put the next plot C (>) YVSTART: (real) Where in the viewport to put the next plot C (>) DELTAX: (real) X size of plot in inches C (>) DELTAY: (real) Y size of plot in inches C (>) MARGIN: (integer) Size of margin in inches C (>) NCOMAX: (integer) Maximum number of coefs C (>) PARAMS: (real array PARAMS(NCOMAX+4,NX,NY)) Final output array C (>) MON: (integer array MON(4,NX,NY)) Monitor array C (>) DATA: (real array DATA(NZ,NKEEP)) Input data array C (>) NKEEP: (integer) Number of spectra kept C (>) NX: (integer) X dimension of the data C (>) NY: (integer) Y dimension of the data C (>) NZ: (integer) Z dimension of the data C (>) XMIN: (real) Lower limit on X to be considered C (>) XMAX: (real) Upper limit on X to be considered C (>) YMIN: (real) Lower limit on Y to be considered C (>) YMAX: (real) Upper limit on Y to be considered C (>) ZMIN: (real) Minimum value of Z C (>) ZMAX: (real) Maximum value of Z C (>) ZAXIS: (real array ZAXIS(NZ)) Z axis data array C (<) NPLOTTED: (integer) Number of spectra plotted C (<) NLAST: (integer) Number of the last spectrum plotted C (>) DLAM: (real) Wavelength increment C (W) WORK1: (real array WORK1(NZ)) Workspace for plotting routines C (W) WORK2: (real array WORK2(NZ)) Workspace for plotting routines C (!) MONITOR: (integer array MONITOR(3,NP*NP)) Plot monitor array C (<) XV1: (real) Min viewport X position of last spectrum plotted C (<) XV2: (real) Max viewport X position of last spectrum plotted C (<) YV1: (real) Min viewport Y position of last spectrum plotted C (<) YV2: (real) Max viewport Y position of last spectrum plotted C (>) SETUP: (logical) TRUE if you are setting up only C (>) REDRAW: (logical) TRUE if most recent panel is to be redrawn. C (>) BACKWARDS: (logical) TRUE if you want the previous panel to be redrawn C C C AUTHOR: Jim Lewis -- RGO C DATE: 28-NOV-1990 C UPDATE: 13-DEC-1990 -- Added ability to go backwards. (JRL) C 18-DEC-1990 -- Changed the order in which spectra are displayed (JRL) C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer np,ncomax,nx,ny,nz,nkeep,nplotted,nlast integer mon(4,nx,ny),monitor(3,np*np) real xvstart,yvstart,deltax,deltay,margin,params(ncomax+4,nx,ny), : zaxis(nz),work1(nz),work2(nz),xmin,xmax,ymin,ymax,zmin,zmax, : dlam,data(nz,nkeep),xv1,xv2,yv1,yv2 logical redraw,backwards,setup C C Functions C integer ich_len C C Local variables C integer ilast,jlast,istart,jstart,tmon,i,j,n,nrem, : ndiv,iadd,dx,xleft,dy,yleft,nxp,ixp,nyp,iyp real fsr_wave,fsr,size,dmin,dmax,xrem(2),yrem(2) character xlab*5,ylab*5,title*12,message*80 C C Common block with FSR C common /free/ fsr C C If this is the first time through, then set up a few variables C if (setup) then C C Set up labels C xlab = ' \gl ' ylab = ' I ' title(1:5) = 'Cell ' C C Set up label sizes C if (np .eq. 1) then size = 1.0 else size = 2.0/float(np) end if C C Set up number of panels depending upon number of spectra included C dx = nint(xmax - xmin) + 1 xleft = mod(dx,np) if (xleft .eq. 0) then nxp = dx/np else nxp = dx/np + 1 end if dy = nint(ymax - ymin) + 1 yleft = mod(dy,np) if (yleft .eq. 0) then nyp = dy/np else nyp = dy/np + 1 end if C C Number of panels plotted so far... C ixp = 0 iyp = 1 go to 500 end if C C Set the character size C call pgsch(size) C C If going backwards, find where the previous panel search C started C If (backwards) then ixp = ixp - 1 if (ixp .le. 0) then iyp = iyp - 1 if (iyp .le. 0) then call dsa_wruser('You CAN''T go backwards!!!') call dsa_wrflush ixp = 0 iyp = 1 go to 500 end if ixp = nxp end if C C If going forwards, then find where in your array was the last plot... C else if (.not. redraw) then ixp = ixp + 1 if (ixp .gt. nxp) then ixp = 1 iyp = iyp + 1 if (iyp .gt. nyp) go to 500 end if end if istart = nint(xmin) + (ixp-1)*np jstart = nint(ymin) + (iyp-1)*np ilast = istart + np - 1 jlast = jstart + np - 1 C C Initialise a few things... C nplotted = 0 tmon = 0 call pgadvance C C FSR in wavelength units... C fsr_wave = fsr*dlam C C Loop through the monitor array and plot each spectrum up individually C yv1 = yvstart do j = jstart,jlast do i = istart,ilast C C If these are outside the X,Y upper limits, then just skip them C (The lower limits are taken care of by definition if istart,jstart) C if ((float(i) .gt. xmax) .or. (float(j) .gt. ymax)) then tmon = tmon + 1 C C Otherwise, continue as normal C else tmon = tmon + 1 n = mon(1,i,j) monitor(1,tmon) = i monitor(2,tmon) = j monitor(3,tmon) = n write(title(6:12),100) i,',',j 100 format(i3,a1,i3) if (tmon .eq. 1) then write(message,101) i,j 101 format('First spectrum on this panel is at',2i4,' ') call dsa_wruser(message(:ich_len(message))) call dsa_wrflush end if C C Find where on the viewport to put this plot C if (np .eq. 1) then xv1 = xvstart yv1 = yvstart else nrem = mod(tmon,np) ndiv = (tmon - 1)/np if (nrem .eq. 1) then xv1 = xvstart yv1 = yvstart + float(ndiv)*(deltay + 2*margin) iadd = 0 else iadd = iadd + 1 xv1 = xvstart + float(iadd)*(deltax + 2*margin) end if end if xv2 = xv1 + deltax yv2 = yv1 + deltay C C Now do the plot. If the spectrum has been left out for one C reason or another (usually it's not between the mask thresholds) C then just put up an empty box C call pgvsize(xv1,xv2,yv1,yv2) if (n .gt. 0) then call minmax(nz,data(1,n),dmin,dmax) else dmin = -1.0 dmax = 1.0 end if call pgwindow(zmin,zmax,dmin,dmax) call pgbox('BCNT',0.0,0,'BCNTS',0.0,0) call pgmtext('B',1.2,0.5,0.5,xlab) call pgmtext('L',1.2,0.5,0.5,ylab) call pgmtext('T',1.2,0.5,0.5,title) if (n .gt. 0) then call ndp_pgbin(zaxis,data(1,n),work1,work2,nz,.false., : .true.) nplotted = nplotted + 1 nlast = n xrem(1) = xv1 xrem(2) = xv2 yrem(1) = yv1 yrem(2) = yv2 end if C C If the fit failed then don't plot the fit. C Otherwise plot a fit based on the fitting function used. C if (n .gt. 0) then if (mon(3,i,j) .gt. 0) then if (mon(2,i,j) .eq. 1) then call gauss_eval(params(1,i,j),ncomax,zaxis,nz, : fsr_wave,work1) else if (mon(2,i,j) .eq. 2) then call cauchy_eval(params(1,i,j),ncomax,zaxis,nz, : fsr_wave,work1) end if C C Now plot the fit as a green line... C call pgsci(3) call pgscr(3,0.0,1.0,0.0) call pgline(nz,zaxis,work1) call pgsci(1) end if end if end if end do end do 500 xv1 = xrem(1) xv2 = xrem(2) yv1 = yrem(1) yv2 = yrem(2) end subroutine neighbour(nx,ny,nz,ncomax,nkeep,fit_type,nco,niter, : boxsize,mon,used,params,data,err,alam0,dlam,errpar,do_it,tried, : resmin,resmax,coefs,errcoefs,picture) C+-------------------------------------------------------------------------------------- C NEIGHBOUR --- Do fits by "nearest neighbour" guessing C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) NX: (integer) X dimension of the data C (>) NY: (integer) Y dimension of the data C (>) NZ: (integer) Z dimension of the data C (>) NCOMAX: (integer) Maximum number of coefs C (>) NKEEP: (integer) Number of spectra kept C (>) FIT_TYPE: (integer) Flag for fitting function C (>) NCO: (integer) Number of coefs C (>) NITER: (integer) Maximum number of iterations for fit. C (>) BOXSIZE: (integer) Binning box size C (>) MON: (integer array MON(4,NX,NY)) Monitor array C (>) USED: (integer array USED(3,NKEEP)) Flags for used spectra C (!) PARAMS: (real array PARAMS(NCOMAX+4,NX,NY)) Final output array C (>) DATA: (real array DATA(NZ,NKEEP)) Input data array C (>) ERR: (real array ERR(NZ)) Input error array C (>) ALAM0: (real) Wavelength zeropoint C (>) DLAM: (real) Wavelength increment C (!) ERRPAR: (real array ERRPAR(NCOMAX+4,NX,NY)) Output array errors C (W) DO_IT: (integer array DO_IT(3,NKEEP)) Mask for pixels with C a neighbour with a good fit. C (W) TRIED: (integer array TRIED(NKEEP)) Flags to tell whether an C attempt has been made here to fit a particular spectrum C (W) RESMIN: (real array RESMIN(NCOMAX)) Minimum values for params C (W) RESMAX: (real array RESMAX(NCOMAX)) Maximum values for params C (W) COEFS: (real array COEFS(NCO)) Temporary storage for coefs C (W) ERRCOEFS: (real array ERRCOEFS(NCO)) Temporary storage for errors C (!) PICTURE: (logical) TRUE if there is a grey scale plot on the screen C C C AUTHOR: Jim Lewis -- RGO C DATE: 28-NOV-1990 C UPDATE: 19-DEC-1990 -- Fitting section taken out (JRL) C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nx,ny,nz,ncomax,nkeep,fit_type,nco,niter,boxsize integer mon(4,nx,ny),used(3,nkeep),do_it(3,nkeep),tried(nkeep) real params(ncomax+4,nx,ny),data(nz,nkeep),err(nz,nkeep), : alam0,dlam,errpar(ncomax+4,nx,ny),coefs(nco),errcoefs(nco), : resmin(ncomax),resmax(ncomax) logical picture C C Functions C integer dsa_typesize integer dyn_element integer ich_len logical par_batch C C Local variables C integer slot,status,ix,iy,bytes_float,bytes_int,zptr,coptr,alptr, : lptr,nfail,nrun,idone,nsweep,i,j,l,nattempt,k,jx,jy,ixrem,iyrem, : address,nspec,totspace,nfail0,ij,ndone,jxst,jxfn,jyst,jyfn, : nditch,ndone0 real dummy,fsr,dum(2),magic_float logical match,ditch,batch character message*80 C C Parameters for DSA and PAR routines... C character range(4)*7,units(4)*7 data range /'rng_amp','rng_pos','rng_sig','rng_bac'/ data units /'counts ','angstrs','angstrs','counts'/ C C Common block with FSR C common /free/ fsr C C Common block for dynamic memory. C include 'DYNAMIC_MEMORY' C C Numeric ranges C include 'NDP_NUMERIC_RANGES' C C Initialise an array which tells whether this routine has already C tried to fit the spectrum... C do i = 1,nkeep tried(i) = 0 end do C C A few variables for workspace allocations C bytes_float = dsa_typesize('float',status) bytes_int = dsa_typesize('int',status) call dsa_get_flag_value('FLOAT',magic_float,status) if (status .ne. 0) go to 500 C C Allocate some workspace for the fitting routines C totspace = (nz + 2*nco*nco)*bytes_float + nco*bytes_int call dsa_get_workspace(totspace,address,slot,status) if (status .ne. 0) then call dsa_wruser : ('Couldn''t allocate workspace in neighbour routine') call dsa_wrflush go to 500 end if zptr = dyn_element(address) coptr = zptr + nz*bytes_float alptr = coptr + nco*nco*bytes_float lptr = alptr + nco*nco*bytes_float C C How many iterations do you want to do? C batch = par_batch() if (.not. batch) call par_cnpar('nsweep') call par_rdval('nsweep',float(min_int),float(max_int),1.0,' ', : dummy) nsweep = nint(dummy) nditch = 0 C C Get limits for coefs C call dsa_wruser : ('Now enter the physical limits for each coef') call dsa_wrflush do i = 1,nco if (.not. batch) call par_cnpar(range(i)) call par_rdary(range(i),min_float,max_float,'a',units(i),2, : 2,dum) resmin(i) = dum(1) resmax(i) = dum(2) end do C C Go through and see how many have already been done. C if (picture) call pgwindow(1.0,float(nx),1.0,float(ny)) ndone = 0 do ij = 1,nkeep i = used(1,ij) j = used(2,ij) nspec = mon(1,i,j) idone = mon(3,i,j) if ((nspec .gt. 0) .and. (idone .gt. 0)) then ndone = ndone + 1 if (picture) call pgpoint(1,float(i),float(j),1) end if end do C C If there aren't any spectra which have been fit, then get out C of here, because you'll have no neighbours! C if (ndone .eq. 0) then call dsa_wruser('No spectra have been fit!!!') call dsa_wrflush go to 500 end if C C Now go through the used array as many times as you've asked for C l = 0 ndone0 = ndone - 1 nattempt = 0 do while ((l .lt. nsweep) .and. (ndone .lt. nkeep) .and. : (ndone0 .ne. ndone)) ndone0 = ndone l = l + 1 C C Go through and find out how many have a neighbour with a valid fit C do ij = 1,nkeep i = used(1,ij) j = used(2,ij) nspec = mon(1,i,j) idone = mon(3,i,j) match = .false. if ((nspec .gt. 0) .and. (idone .le. 0) .and. : (tried(ij) .eq. 0)) then ix = i iy = j jxst = max(1,ix-1) jxfn = min(nx,ix+1) jyst = max(1,iy-1) jyfn = min(ny,iy+1) jy = jyst - 1 do while ((.not. match) .and. (jy .lt. jyfn)) jy = jy + 1 jx = jxst - 1 do while ((.not. match) .and. (jx .lt. jxfn)) jx = jx + 1 if (mon(3,jx,jy) .gt. 0) then match = .true. ixrem = jx iyrem = jy end if end do end do end if C C If this pixel does have a valid neighbour, then set up some flags C if (match) then do_it(1,ij) = 1 do_it(2,ij) = ixrem do_it(3,ij) = iyrem else do_it(1,ij) = 0 end if end do C C Now that you've gone through and marked all those with a valid C neighbour, try and do the fits... C do ij = 1,nkeep if (do_it(1,ij) .eq. 1) then ix = used(1,ij) iy = used(2,ij) nspec = mon(1,ix,iy) nattempt = nattempt + 1 ndone = ndone + 1 ixrem = do_it(2,ij) iyrem = do_it(3,ij) coefs(1) = params(1,ixrem,iyrem) coefs(2) = (params(2,ixrem,iyrem) - alam0)/dlam coefs(3) = params(3,ixrem,iyrem)/abs(dlam) coefs(4) = params(4,ixrem,iyrem) C C Do the fit. C nfail0 = nfail call fit_spec(nco,ncomax,dynamic_mem(lptr),nfail, : nrun,nx,ny,nz,0,nkeep,nspec, : data(1,nspec),err(1,nspec),used,dynamic_mem(zptr), : fsr,coefs,errcoefs,niter,dynamic_mem(coptr), : dynamic_mem(alptr),fit_type,alam0,dlam,params) if (picture) call pgpoint(1,float(ix),float(iy),1) C C If the fit failed, then flag it... C if (nfail .ne. nfail0) then do k = 1,ncomax params(k,ix,iy) = magic_float errpar(k,ix,iy) = magic_float end do mon(2,ix,iy) = fit_type mon(3,ix,iy) = -3 tried(nspec) = 1 C C Otherwise, add it into the parameters list C else ditch = .false. do k = 1,nco if ((coefs(k) .lt. resmin(k)) .or. : (coefs(k) .gt. resmax(k))) ditch = .true. end do if (ditch) then do k = 1,ncomax params(k,ix,iy) = magic_float errpar(k,ix,iy) = magic_float end do mon(2,ix,iy) = fit_type mon(3,ix,iy) = -1 nditch = nditch + 1 tried(nspec) = 1 else do k = 1,nco params(k,ix,iy) = coefs(k) errpar(k,ix,iy) = errcoefs(k) end do mon(2,ix,iy) = fit_type mon(3,ix,iy) = 1 mon(4,ix,iy) = boxsize tried(nspec) = 1 end if end if end if end do write(message,100) l,nattempt,ndone,nkeep 100 format(i4,' sweeps',i7,' spectra attempted',i7,' done out of', : i7) call dsa_wruser(message(:ich_len(message))) call dsa_wrflush end do write(message,101) nfail 101 format(i7,' spectra failed') call dsa_wruser(message(:ich_len(message))) call dsa_wrflush write(message,102) nrun 102 format(i7,' spectra had runaway solutions') call dsa_wruser(message(:ich_len(message))) call dsa_wrflush write(message,105) nditch 105 format(i7,' spectra had unphysical solutions') call dsa_wruser(message(:ich_len(message))) call dsa_wrflush C C Tidy and exit C 500 if (slot .ne. 0) call dsa_free_workspace(slot,status) slot = 0 end subroutine info(fit_type,iwt,how_guess,boxsize,niter, : ntriv_min,dchi_max) C+-------------------------------------------------------------------------------------- C INFO --- Tell the user some useful information about his/her/its setup C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) FIT_TYPE: (integer) Flag for the fitting function to be used C (>) IWT: (integer) Tells which weighting scheme is being used C (>) HOW_GUESS: (integer) Flag for the guessing algorithm used C (>) BOXSIZE: (integer) Size of the smoothing box used. C (>) NITER: (integer) Maximum number of iterations to be done C (>) NTRIV_MIN: (integer) Minimum number of trivial changes for convergence C (>) DCHI_MAX: (real) Definition of trivial fractional change in chi-squared C C C AUTHOR: Jim Lewis -- RGO C DATE: 18-DEC-1990 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer fit_type,iwt,how_guess,boxsize,niter,ntriv_min real dchi_max C C Functions C integer ich_len C C Local variables C real yrange(2) character message*80 logical scale_y C C Common block with Y scaling info for line plotting routines C common /scaling/ scale_y,yrange C C Write out fit type... C call dsa_wruser(' ') call dsa_wrflush if (fit_type .eq. 1) then message = 'gaussian' else message = 'cauchy function' end if call dsa_wruser('The fitting function is a '// : message(:ich_len(message))) call dsa_wrflush C C Write out weighting scheme C if (iwt .eq. 1) then message = 'sqrt(n) weights' else message = 'equal weights' end if call dsa_wruser('The points are being weighted by '// : message(:ich_len(message))) call dsa_wrflush C C Write out guessing algorithm C if (how_guess .eq. 1) then message = 'moments analysis' else if (how_guess .eq. 2) then message = 'previous results' else if (how_guess .eq. 3) then message = 'hand (slow version)' else if (how_guess .eq. 4) then message = 'hand (fast version)' end if call dsa_wruser('The guessing is done by '// : message(:ich_len(message))) call dsa_wrflush C C Write out the smoothing box size C write(message,100) boxsize 100 format('The smoothing box size is:',i4,' pixels') call dsa_wruser(message(:ich_len(message))) call dsa_wrflush C C Write out maximum number of iterations C write(message,101) niter 101 format('The maximum number of iterations to be done is:',i4) call dsa_wruser(message(:ich_len(message))) call dsa_wrflush C C Write out the number of trivial changes for convergence C write(message,102) ntriv_min 102 format('The number of trivial changes for convergence is:', : i4) call dsa_wruser(message(:ich_len(message))) call dsa_wrflush C C Write out the definition of a trivial fractional change in chi-squared C write(message,103) dchi_max 103 format('Trivial fractional change in chi-squared defined as:', : f8.4) call dsa_wruser(message(:ich_len(message))) call dsa_wrflush C C Finally, what scale are the line plots on? C if (scale_y) then write(message,104) yrange(1),yrange(2) 104 format('Y axis on all plots will have',2f8.3, : ' as Min and Max') call dsa_wruser(message(:ich_len(message))) call dsa_wrflush else call dsa_wruser('All plots will be scaled individually') call dsa_wrflush end if call dsa_wruser(' ') call dsa_wrflush end subroutine cursor_point(nx,ny,xaxis,yaxis,xid,yid,quit) C+-------------------------------------------------------------------------------------- C CURSOR_POINT --- Find a single point with the cursor C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) NX: (integer) X dimension of data C (>) NY: (integer) Y dimension of data C (>) XAXIS: (real array XAXIS(NX)) X axis data C (>) YAXIS: (real array YAXIS(NY)) Y axis data C (<) XID: (real) Output X coordinate (pixel units) C (<) YID: (real) Output Y coordinate (pixel units) C (<) QUIT: (logical) TRUE if user quit rather than indicating a point C C C AUTHOR: Jim Lewis -- RGO C DATE: 19-DEC-1990 C UPDATE: C C+-------------------------------------------------------------------------------------- implicit none C C Parameters C integer nx,ny real xaxis(nx),yaxis(ny),xid,yid logical quit C C Local variables C integer stapix(2),endpix(2),nid real start(2),end(2),ximv(2),yimv(2),square,xlast,ylast C C Set up axes C stapix(1) = 1 stapix(2) = 1 endpix(1) = nx endpix(2) = ny start(1) = xaxis(1) start(2) = yaxis(1) end(1) = xaxis(nx) end(2) = yaxis(ny) C C Choose the points C call ndp_image_viewport(stapix,endpix,1.0,'L',ximv,yimv, : square) quit = .false. call ndp_image_cursor(xaxis,yaxis,nx,ny,start,end, : stapix,endpix,1,'SDP',ximv,yimv,xid,yid,nid,xlast, : ylast,quit) end subroutine cleanvel(nx,ny,ncomax,mon,params,xaxis,yaxis,fsr,alam0, : dlam,picture,ref_name,mask,work1,work2,work3,work4) C+-------------------------------------------------------------------------------------- C CLEANVEL --- Remove jumps of 1 or more FSRs from velocity map. C C Parameters (">" Input, "<" Output, "W" Workspace, "!" Changed) C C (>) NX: (integer) X dimension of the data C (>) NY: (integer) Y dimension of the data C (>) NCOMAX: (integer) Maximum number of coefs C (>) MON: (integer array MON(4,NX,NY)) Monitor array C (>) PARAMS: (real array PARAMS(NCOMAX+4,NX,NY)) Final output array C (>) XAXIS: (real array XAXIS(NX)) X axis data array C (>) YAXIS: (real array YAXIS(NY)) Y axis data array C (>) FSR: (real) Free spectral range C (>) ALAM0: (real) Wavelength zeropoint C (>) DLAM: (real) Wavelength increment C (>) PICTURE: (logical) TRUE if a 2-d plot is on the screen C (>) REF_NAME: (character) Reference name of input cube C (W) MASK: (byte array MASK(NX,NY)) Workspace C (W) WORK1: (real array WORK1(NX,NY)) Workspace C (W) WORK2: (integer array WORK2(NX,NY)) Workspace C (W) WORK3: (real array WORK3(NX,NY)) Workspace C (W) WORK4: (integer array WORK4(NX,NY)) Workspace C C C AUTHOR: Jim Lewis -- RGO C DATE: 25-FEB-1991 C UPDATE: 29-JAN-1992 -- Modified to use new NDP 2d plot routines (JRL) C C+-------------------------------------------------------------------------------------- implicit none C C Paramters C integer nx,ny,mon(4,nx,ny),ncomax,work2(nx,ny),work4(nx,ny) real params(ncomax+4,nx,ny),work1(nx,ny),xaxis(nx),yaxis(nx),fsr, : alam0,dlam,work3(nx,ny) byte mask(nx,ny) logical picture character ref_name*(*) C C Functions C integer dsa_typesize integer dyn_element integer ich_len logical par_batch logical par_qnum logical par_quest C C Local variables C integer i,j,idone,status,ignore2,stapix(2),endpix(2),jx,jy, : jxst,jxfn,jyst,jyfn,ixrem,iyrem,numchanged,isweep, : nsweep,vptr,vslot,address,odims(2),bytes_float,bytes_int, : ix,iy,nord,ndim,nelm real lam_cen,c,lamda,vel,low,high,start(2),end(2),dv, : xx,yy,tol,pct_tol,dummy,vmin,vmax,fsr_vel, : lownew,highnew,magic_float character message*80,oldtable*20,table*20,vel_field*80,type*20 logical badpix,continue,changed,quit,match,accept,writeout,ignore, : looked,already,batch C C Common block for dynamic memory C include 'DYNAMIC_MEMORY' C C Numeric ranges C include 'NDP_NUMERIC_RANGES' C C Initialize a few things C status = 0 changed = .false. badpix = .false. vmin = max_float vmax = min_float bytes_float = dsa_typesize('float',status) bytes_int = dsa_typesize('int',status) call dsa_get_flag_value('FLOAT',magic_float,status) if (status .ne. 0) go to 500 C C Set FSR in km/sec C c = 3.0e5 fsr_vel = c*fsr*dlam/alam0 C C Find out if you already have a velocity map C already = par_quest('Do you already have a velocity map?',.true.) C C If you already have a velocity map, then open the file for it and C move the information in it to the first workspace C batch = par_batch() if (already) then if (.not. batch) call par_cnpar('vel_field') call par_rdchar('vel_field/nocheck',' ',vel_field) call dsa_named_input('vel_field',vel_field,status) call ndp_get_image_info('vel_field',.true.,.false., : type,badpix,status) call dsa_data_size('vel_field',2,ndim,odims,nelm,status) if (status .ne. 0) go to 500 if (ndim .ne. 2) then call dsa_wruser('Plot image must be 2-d') call dsa_wrflush call dsa_close_structure('vel_field',status) go to 500 end if if ((odims(1) .ne. nx) .and. (odims(2) .ne. ny)) then call dsa_wruser : ('Velocity field image must be same size as data cube') call dsa_wrflush call dsa_close_structure('vel_field',status) go to 500 end if if (badpix) call dsa_use_flagged_values('vel_field',status) call dsa_map_data('vel_field','read','float',address,vslot, : status) if (status .ne. 0) go to 500 vptr = dyn_element(address) call gen_move(nx*ny*bytes_float,dynamic_mem(vptr),work1) if (badpix) : call dsa_post_process_flagged_values('vel_field',status) call dsa_unmap(vslot,status) call dsa_close_structure('vel_field',status) if (status .ne. 0) go to 500 C C Find min and max of field. C do j = 1,ny do i = 1,nx if (mon(3,i,j) .gt. 0) then vmin = min(vmin,work1(i,j)) vmax = max(vmax,work1(i,j)) end if end do end do C C Otherwise, calculate a velocity field from the parameters C else C C First get central wavelength C call dsa_wruser('OK, I''ll generate one...') call dsa_wrflush if (.not. batch) call par_cnpar('lam_cen') call par_rdval('lam_cen',min_float,max_float,1.0,'angstroms', : lam_cen) C C Now calculate velocity field and put it in the velocity workspace C do j = 1,ny do i = 1,nx idone = mon(3,i,j) if (idone .gt. 0) then lamda = params(2,i,j) vel = c*(lamda - lam_cen)/lam_cen work1(i,j) = vel vmin = min(vmin,vel) vmax = max(vmax,vel) else work1(i,j) = magic_float badpix = .true. end if end do end do end if C C Tell the user the ranges of the data and ask for new ones C write(message,100) vmin,vmax 100 format('Range of velocity is:',2f12.2,' km/sec') call dsa_wruser(message(:ich_len(message))) call dsa_wrflush ignore = par_qnum('Low value for plot',min_float, : max_float,vmin,.true.,' ',low) ignore = par_qnum('High value for plot',min_float, : max_float,vmax,.true.,' ',high) call var_getchr('table',0,0,oldtable,ignore2) call par_qstr('Name of lookup table',oldtable,.false., : .false.,table) call var_setchr('table',0,0,table,ignore2) C C Set up colour index map C call ndp_image_index(nx*ny,low,high,work1,badpix,work2,status) if (status .ne. 0) go to 450 C C Set up axes C stapix(1) = 1 stapix(2) = 1 endpix(1) = nx endpix(2) = ny start(1) = xaxis(1) start(2) = yaxis(1) end(1) = xaxis(nx) end(2) = yaxis(ny) C C Now display the first 2-d plot of the velocity field C call disp(work2,nx,ny,stapix,endpix,start,end,high,low, : picture,table,status) if (status .ne. 0) go to 500 C C Now loop so long as you really want C continue = par_quest('Clean this map?',.true.) do while (continue) C C Select pixel to serve as fiducial pixel C call dsa_wruser : ('Indicate the fiducial pixel with the cursor') call dsa_wrflush call cursor_point(nx,ny,xaxis,yaxis,xx,yy,quit) if (quit) go to 450 ix = nint(xx) iy = nint(yy) if (mon(3,ix,iy) .le. 0) then call dsa_wruser('This pixel doesn''t have a valid fit') call dsa_wrflush go to 450 end if C C Zero mask image C call gen_fill(nx*ny*bytes_int,0,mask) mask(ix,iy) = 1 C C Now set mask image for pixels with no fit to 3 C do j = 1,ny do i = 1,nx if (mon(3,i,j) .le. 0) mask(i,j) = 3 end do end do C C Get tolerance and number of sweeps C if (.not. batch) call par_cnpar('tol') call par_rdval('tol',min_float,max_float,1.0,'percentage', : pct_tol) tol = 0.01*pct_tol*fsr_vel if (.not. batch) call par_cnpar('nsweep') call par_rdval('nsweep',float(min_int),float(max_int),1.0, : ' ',dummy) nsweep = nint(dummy) C C Copy the old map... C call gen_move(nx*ny*bytes_float,work1,work3) C C Now do the correction C isweep = 0 looked = .true. numchanged = 0 do while ((looked) .and. (isweep .lt. nsweep)) isweep = isweep + 1 C C Loop through map and check for pixels with a good fit C looked = .false. do j = 1,ny do i = 1,nx if (mask(i,j) .eq. 0) then jxst = max(1,i-1) jxfn = min(nx,i+1) jyst = max(1,j-1) jyfn = min(ny,j+1) C C Loop through cell to find match C jy = jyst - 1 match = .false. do while ((.not. match) .and. (jy .lt. jyfn)) jy = jy + 1 jx = jxst - 1 do while ((.not. match) .and. (jx .lt. jxfn)) jx = jx + 1 if (mask(jx,jy) .eq. 1) then match = .true. ixrem = jx iyrem = jy end if end do end do C C Now check on the order number of the velocity if there C has been a match C if (match) then looked = .true. dv = work3(ixrem,iyrem) - work1(i,j) nord = nint(dv/fsr_vel) dv = dv - nord*fsr_vel if ((abs(dv) .le. tol) .and. (nord .ne. 0)) then vel = work1(i,j) + nord*fsr_vel work3(i,j) = vel numchanged = numchanged + 1 end if mask(i,j) = 2 end if end if end do end do C C Reset all changed and otherwise good velocities so that C they can be used as fiducial spectra in the next sweep C do j = 1,ny do i = 1,nx if (mask(i,j) .eq. 2) mask(i,j) = 1 end do end do write(message,101) isweep,numchanged 101 format(i4,' sweeps',i7,' velocities switched') call dsa_wruser(message(:ich_len(message))) call dsa_wrflush end do C C Find min and max of new velocity array C vmin = max_float vmax = min_float do j = 1,ny do i = 1,nx if (mon(3,i,j) .gt. 0) then vmin = min(work3(i,j),vmin) vmax = max(work3(i,j),vmax) end if end do end do C C Write out new min and max and get min and max for new 2-d plot C write(message,102) vmin,vmax 102 format('New range of velocity is:',2f12.2,' km/sec') call dsa_wruser(message(:ich_len(message))) call dsa_wrflush ignore = par_qnum('Low value for plot',min_float, : max_float,vmin,.true.,' ',lownew) ignore = par_qnum('High value for plot',min_float, : max_float,vmax,.true.,' ',highnew) C C Set up new colour index map C call ndp_image_index(nx*ny,lownew,highnew,work3,badpix, : work4,status) if (status .ne. 0) go to 450 C C Now display the revised 2-d plot of the velocity field C call disp(work4,nx,ny,stapix,endpix,start,end,highnew,lownew, : picture,table,status) if (status .ne. 0) then status = 0 go to 450 end if C C See if the user wants to accept these changes C accept = par_quest('Accept these changes?',.true.) C C If the user accepts the changes, then move the new results into C the old results workspace and continue. Otherwise plot up old C results and continue C if (accept) then call gen_move(nx*ny*bytes_float,work3,work1) high = highnew low = lownew changed = .true. else call dsa_wruser('Plotting up old field...') call dsa_wrflush call disp(work2,nx,ny,stapix,endpix,start,end,high,low, : picture,table,status) end if C C See if the user wants to do any more... C continue = par_quest('Try to clean this map again?',.true.) end do C C If any changes have occured in the velocity field, then write it out C if you so desire. C 450 if (changed) then call dsa_wruser('Some changes have been made to the velocity') call dsa_wruser(' field') call dsa_wrflush writeout = par_quest('Do you want to write this to a file?', : .true.) if (writeout) then if (.not. batch) call par_cnpar('vel_field') call par_rdchar('vel_field/nocheck',' ',vel_field) call dsa_named_output('vel_field',vel_field,' ',1,1, : status) odims(1) = nx odims(2) = ny call dsa_simple_output('vel_field','d,a1,a2','float',2, : odims,status) call dsa_reshape_axis('vel_field',1,ref_name,2,1,nx, : status) call dsa_reshape_axis('vel_field',2,ref_name,3,1,ny, : status) call dsa_set_flagged_values('vel_field',.true.,status) call dsa_use_flagged_values('vel_field',status) call dsa_map_data('vel_field','write','float',address, : vslot,status) if (status .ne. 0) go to 500 vptr = dyn_element(address) call gen_move(nx*ny*bytes_float,work1,dynamic_mem(vptr)) call dsa_set_range('vel_field',vmin,vmax,status) call dsa_post_process_flagged_values('vel_field',status) call dsa_unmap(vslot,status) call dsa_close_structure('vel_field',status) end if end if 500 end