Merge prs

build/.depends

@@ -74,7 +74,7 @@ libread.o : libread.f90 libutils.o globaldata.o globaldata.o globaldata.o asciiu
 libutils.o : libutils.F90 libinclude.f90 
 lightcurve.o : lightcurve.f90 write_fits.o system_utils.o rotate.o options_xsecrotate.o blackbody.o globaldata.o write_pixmap.o physcon.o units.o options_render.o options_particleplots.o globaldata.o interpolation.o globaldata.o kernels.o interpolate3D_opacity.o lightcurve_utils.o limits.o labels.o globaldata.o 
 lightcurve_utils.o : lightcurve_utils.f90 utils_vectors.o physcon.o globaldata.o 
-limits.o : limits.f90 asciiutils.o globaldata.o globaldata.o geometry.o labels.o globaldata.o 
+limits.o : limits.f90 system_utils.o asciiutils.o globaldata.o globaldata.o geometry.o labels.o globaldata.o 
 map_columns.o : map_columns.f90 prompting.o promptlist.o labels.o asciiutils.o 
 menu.o : menu.f90 write_pixmap.o system_utils.o geometry.o globaldata.o globaldata.o timestepping.o geomutils.o get_data.o defaults.o transform.o prompting.o globaldata.o units.o options_xsecrotate.o options_vecplot.o options_render.o options_page.o options_particleplots.o options_limits.o globaldata.o options_data.o limits.o labels.o globaldata.o 
 moments.o : moments.f90 timing.o interpolate1D.o 
@@ -150,7 +150,7 @@ setpage.o : setpage.f90 transform.o asciiutils.o ${PLOTLIB:.f90=.o}
 shapes.o : shapes.f90 legends.o interactive_help.o parsetext.o asciiutils.o transform.o ${PLOTLIB:.f90=.o} exact_function.o prompting.o globaldata.o promptlist.o 
 sort.o : sort.f90 
 sph_moments.o : sph_moments.f90 kernels.o moments.o system_utils.o write_fits.o 
-splash.o : splash.f90 globaldata.o interactive.o labels.o globaldata.o exact.o set_options_from_dataread.o read_data.o colours.o options_xsecrotate.o options_render.o options_particleplots.o options_page.o timestepping.o analysis.o write_griddata.o write_sphdata.o convert.o write_pixmap.o asciiutils.o system_utils.o system_f2003.o globaldata.o interpolate3D_projection.o allocate.o menu.o kernels.o limits.o initialise.o defaults.o geomutils.o get_data.o globaldata.o 
+splash.o : splash.f90 globaldata.o interactive.o labels.o globaldata.o exact.o set_options_from_dataread.o read_data.o colours.o options_xsecrotate.o options_render.o options_particleplots.o options_limits.o options_page.o timestepping.o analysis.o write_griddata.o write_sphdata.o convert.o write_pixmap.o asciiutils.o system_utils.o system_f2003.o globaldata.o interpolate3D_projection.o allocate.o menu.o kernels.o limits.o initialise.o defaults.o geomutils.o get_data.o globaldata.o 
 system_f2003.o : system_f2003.f90 
 system_utils.o : system_utils.f90 system_f2003.o asciiutils.o 
 timestepping.o : timestepping.f90 globaldata.o options_particleplots.o get_data.o ${PLOTLIB:.f90=.o} plotstep.o options_page.o globaldata.o globaldata.o globaldata.o 

docs/getting-started.rst

@@ -442,10 +442,27 @@ Command line flags (or environment variables) that affect all data reads are:
 |                      |                       | :ref:`sec:menu-r` is also turned on, otherwise  |
 |                      |                       | slow rendering can result.                      |
 +----------------------+-----------------------+-------------------------------------------------+
-| ---xmin=0.1          | SPLASH_MARGIN_XMIN    | can be used to manually adjust the left page    |
-| ---xmax=0.1          | SPLASH_MARGIN_XMAX    | page margin (set to fraction of viewport,       |
-| ---ymin=0.1          | SPLASH_MARGIN_YMIN    | negative values are allowed).                   |
-| ---ymax=0.1          | SPLASH_MARGIN_YMAX    |                                                 |
+| ---xmin=0.0          |                       | override plot limits for x,y,z coordinate       |
+| ---xmax=1.0          |                       | columns after splash.limits is read.            |
+| ---ymin=0.0          |                       | Mutually exclusive with --limits and --lim.     |
+| ---ymax=1.0          |                       |                                                 |
+| ---zmin=0.0          |                       |                                                 |
+| ---zmax=1.0          |                       |                                                 |
++----------------------+-----------------------+-------------------------------------------------+
+| ---limits=-1.0,1.0   |                       | set plot limits:                                |
+| ---limits=[-1.0,1.0] |                       | - 2 values: same range on all coordinate axes   |
+|                      |                       | - 2*ndim values: one min,max pair per axis      |
+|                      |                       | - other counts: one min,max pair per column     |
+|                      |                       | Mutually exclusive with --xmin and --lim.       |
++----------------------+-----------------------+-------------------------------------------------+
+| ---lim=1.0           |                       | centred box [-lim,lim] on all coordinate axes   |
+|                      |                       | (same as menu "centred cube"). Exclusive with   |
+|                      |                       | --xmin and --limits.                            |
++----------------------+-----------------------+-------------------------------------------------+
+| ---xminmargin=0.1    | SPLASH_XMINMARGIN     | manually adjust the page margins (fractions     |
+| ---xmaxmargin=0.1    | SPLASH_XMAXMARGIN     | of viewport; negative values allowed).          |
+| ---yminmargin=0.1    | SPLASH_YMINMARGIN     |                                                 |
+| ---ymaxmargin=0.1    | SPLASH_YMAXMARGIN     |                                                 |
 +----------------------+-----------------------+-------------------------------------------------+
 
 .. _sec:splash:

docs/lightcurve.rst

@@ -74,11 +74,14 @@ where the :math:`w_{\mathrm{col}, ij}` is the dimensionless column kernel.
 
 The spectra were then summed across the frequencies and pixels to produce the bolometric luminosity, listed in [02  Luminosity] in the outputted .out file.
 
+The bounds of the ray-traced region are taken from the plot limits for the coordinate columns (from the particle extent, or from ``splash.limits`` if present). These can be overridden with the global flags ``--xmin``, ``--xmax``, ``--ymin``, ``--ymax``, ``--limits=...`` or ``--lim=...`` (see *Options affecting all data reads* in Getting started).
 
 When ``--temperature`` is set, temperature :math:`T` is computed at **read time** from density :math:`\rho` and
 specific internal energy :math:`u` using *get_temp_from_u()* in ``lightcurve_utils.f90``, assuming gas and
 radiation pressure in local thermodynamic equilibrium (see :ref:`sec:lckappa`).
 
+The temperature for each SPH particle can be altered with a scaling factor (the spectral hardening factor) and an optional constant offset. The scaling factor is set with ``--fcol`` (or ``--f_col``) and the offset with ``--ocol``. The scaling is applied first: :math:`T \rightarrow f_{\rm col} T + o_{\rm col}`.
+
 Opacity :math:`\kappa` used in the ray trace is described in :ref:`sec:lckappa`.
 
 

src/convert.f90

@@ -88,6 +88,7 @@ subroutine convert_all(outformat,igotfilenames,useall)
        call get_data(ifile,igotfilenames,firsttime=.true.)
        !
        ! read plot limits from file (overrides get_data limits settings)
+       ! and apply any command-line limit overrides (--xmin, --limits, --lim, ...)
        !
        call read_limits(trim(limitsfile),ierr)
        !

src/lightcurve.f90

@@ -83,7 +83,7 @@ subroutine get_lightcurve(ncolumns,dat,npartoftype,masstype,itype,ndim,ntypes,&
  real, dimension(:),   allocatable :: freq,spectrum,bb_spectrum
  real, dimension(:,:), allocatable :: img,taupix,flux_nu,v_on_c,badpix
  real, dimension(:,:,:), allocatable :: img_nu,img_tmp
- real :: zobs,dzobs,dx,dy,area,freqmin,freqmax,lam_max,freq_max,bb_scale,opacity_factor,f_col
+ real :: zobs,dzobs,dx,dy,area,freqmin,freqmax,lam_max,freq_max,bb_scale,opacity_factor,f_col,o_col
  real :: betaz,lorentz,doppler_factor,doppler_factor_max,tempi,badarea
  real :: rstar,lstar
  logical :: relativistic,nofits
@@ -127,8 +127,18 @@ subroutine get_lightcurve(ncolumns,dat,npartoftype,masstype,itype,ndim,ntypes,&
  if (f_col < 0.0) f_col = renvironment('f_col', 1.0)
  print "(/,a,f5.2,/,a,/)",' SPECTRAL HARDENING FACTOR f_col = ',f_col, &
                           ' (use --fcol=1.7 for radiation-pressure dominated flows; Shimura & Takahara 1995)'
+ o_col = renvironment('ocol',0.)
+ if (abs(o_col) > tiny(0.)) print "(a,1pg10.3)",' temperature offset o_col = ',o_col
+
+
+ !--render region from plot limits (splash.limits and/or --xmin/--limits/--lim)
  xmin(1:ndim) = lim(ix(1:ndim),1)
  xmax(1:ndim) = lim(ix(1:ndim),2)
+ if (xmax(1) <= xmin(1) .or. xmax(2) <= xmin(2)) then
+    print "(a)",' ERROR: invalid render bounds (check .limits file or command line flags)'
+    ierr = 3
+    return
+ endif
  !
  !--set number of particles to use in the interpolation routines
  !  and allocate memory for weights
@@ -137,7 +147,7 @@ subroutine get_lightcurve(ncolumns,dat,npartoftype,masstype,itype,ndim,ntypes,&
  allocate(weight(n),x(n),y(n),z(n),flux(n),opacity(n),h(n),stat=ierr)
  if (ierr /= 0) then
     print*,' ERROR allocating memory for interpolation weights, aborting...'
-    ierr = 3
+    ierr = 4
     return
  endif
  x(1:n) = dat(1:n,ix(1))
@@ -198,11 +208,8 @@ subroutine get_lightcurve(ncolumns,dat,npartoftype,masstype,itype,ndim,ntypes,&
     print "(a,1pg10.2,a)",' WARNING: using fixed opacity kappa = ',maxval(opacity),' cm^2/g for lightcurve'
  endif
  !
- ! specify source function for each particle
+ ! specify source function for each particle (can be frequency-dependent)
  !
- flux = steboltz*dat(1:n,itemp)**4 ! grey version
-
- ! frequency-dependent version
  nfreq = 128
  freqmin = 1e8
  freqmax = 1e22
@@ -214,7 +221,7 @@ subroutine get_lightcurve(ncolumns,dat,npartoftype,masstype,itype,ndim,ntypes,&
  doppler_factor_max = 0.
  !$omp parallel do default(none) &
  !$omp shared(n,nfreq,freq,flux,flux_nu,dat,h,weight,radkernel) &
- !$omp shared(opacity,relativistic,v_on_c,itemp,f_col) &
+ !$omp shared(opacity,relativistic,v_on_c,itemp,f_col,o_col) &
  !$omp private(i,betaz,lorentz,tempi,rstar,lstar) &
  !$omp firstprivate(opacity_factor,doppler_factor) &
  !$omp reduction(max:doppler_factor_max)
@@ -239,12 +246,19 @@ subroutine get_lightcurve(ncolumns,dat,npartoftype,masstype,itype,ndim,ntypes,&
        lstar = 4.*pi*rstar**2*steboltz*tempi**4
        print "(a,2(es10.3,a),/)",' Luminosity of sink = ',lstar,' erg/s = ',lstar/Lsun,' L_sun'
     else
-       tempi = dat(i,itemp)*f_col
+       tempi = dat(i,itemp)*f_col + o_col
     endif
     !call get_opacity_nongrey(nfreq,freq,dat(i,temp),dat(i,rho),opacity_nu(:,i))
-    flux_nu(:,i) = B_nu(tempi,freq*doppler_factor)
+    if (tempi > 0.) then
+       flux_nu(:,i) = B_nu(tempi,freq*doppler_factor)
+       flux(i) = steboltz*tempi**4 ! grey version
+    else
+       flux_nu(:,i) = 0.
+       flux(i) = 0.
+    endif
  enddo
  !$omp end parallel do
+
  if (relativistic) print*,' max relativistic correction=',doppler_factor_max
 
  if (allocated(img_nu)) deallocate(img_nu)

src/limits.f90

@@ -31,7 +31,8 @@ module limits
  use params
  implicit none
  real, dimension(maxplot,2) :: lim,range,lim2
- private :: warn_minmax
+ private :: warn_minmax,check_coord_column,set_coord_limits,&
+            read_limits_file,override_limits_from_flags
 
  public
 
@@ -154,10 +155,23 @@ subroutine write_limits(limitsfile,iverbose)
 
 end subroutine write_limits
 
+!----------------------------------------------------------
+! read plot limits from file, then apply command-line overrides
+!----------------------------------------------------------
+subroutine read_limits(limitsfile,ierr,overridden)
+ character(len=*), intent(in) :: limitsfile
+ integer,         intent(out) :: ierr
+ logical, intent(out), optional :: overridden
+
+ call read_limits_file(limitsfile,ierr)
+ call override_limits_from_flags(overridden)
+
+end subroutine read_limits
+
 !----------------------------------------------------------
 ! read plot limits for all columns from a file
 !----------------------------------------------------------
-subroutine read_limits(limitsfile,ierr)
+subroutine read_limits_file(limitsfile,ierr)
  use labels,        only:label
  use settings_data, only:numplot,ncolumns,ncalc,iverbose
  use asciiutils,    only:ncolumnsline
@@ -227,7 +241,7 @@ subroutine read_limits(limitsfile,ierr)
  endif
  close(unit=54)
 
-end subroutine read_limits
+end subroutine read_limits_file
 
 !----------------------------------------------------------
 ! get a subset of the particles by enforcing range restrictions
@@ -468,4 +482,140 @@ logical function limits_are_equal(n,iplotx,iploty)
 
 end function limits_are_equal
 
+!----------------------------------------------------------
+! check that coordinate axis idim maps to a valid plot column
+!----------------------------------------------------------
+subroutine check_coord_column(idim,flag)
+ use labels,        only:ix
+ use settings_data, only:ndim,numplot
+ integer, intent(in) :: idim
+ character(len=*), intent(in) :: flag
+
+ if (idim < 1 .or. idim > ndim .or. ix(idim) <= 0 .or. ix(idim) > numplot) then
+    print "(a)",' ERROR: '//trim(flag)//' given but coordinate column is not set'
+    stop
+ endif
+
+end subroutine check_coord_column
+
+!----------------------------------------------------------
+! set lim for coordinate axis idim after validating ix(idim)
+!----------------------------------------------------------
+subroutine set_coord_limits(idim,xlo,xhi,flag)
+ use labels, only:ix
+ integer, intent(in) :: idim
+ real, intent(in) :: xlo,xhi
+ character(len=*), intent(in) :: flag
+
+ call check_coord_column(idim,flag)
+ lim(ix(idim),1:2) = (/xlo,xhi/)
+
+end subroutine set_coord_limits
+
+!----------------------------------------------------------
+! apply --xmin/--limits/--lim overrides (mutually exclusive)
+!----------------------------------------------------------
+subroutine override_limits_from_flags(overridden)
+ use labels,        only:ix
+ use geometry,      only:labelcoord
+ use settings_data, only:ndim,numplot
+ use system_utils,  only:get_command_flag,get_command_option,rflaglist
+ integer, parameter :: maxvals = 2*maxplot
+ real :: vals(maxvals),xmin,xmax,halfwidth
+ logical, intent(out), optional :: overridden
+ logical :: have_axis,have_limits,have_lim,have_min(3),have_max(3)
+ integer :: i,nvals,nfam,icol
+
+ do i=1,3
+    have_min(i) = get_command_flag(trim(labelcoord(i,1))//'min')
+    have_max(i) = get_command_flag(trim(labelcoord(i,1))//'max')
+ enddo
+ have_axis = any(have_min .or. have_max)
+ have_limits = get_command_flag('limits')
+ have_lim    = get_command_flag('lim')
+
+ nfam = count((/have_axis,have_limits,have_lim/))
+ if (present(overridden)) overridden = (nfam > 0)
+ if (nfam > 1) then
+    print "(a)",' ERROR: use only one of --xmin/--xmax/... or --limits=... or --lim=...'
+    stop
+ elseif (nfam == 0) then
+    return
+ endif
+
+ if (have_lim) then
+    !
+    ! set x, y and z limits from --lim=halfwidth flag
+    !
+    halfwidth = get_command_option('lim',default=-1.)
+    if (halfwidth <= 0.) then
+       print "(a)",' ERROR: --lim requires a positive halfwidth, e.g. --lim=1.0'
+       stop
+    elseif (ndim < 1) then
+       print "(a)",' ERROR: --lim given but coordinate columns are not set'
+       stop
+    endif
+    do i=1,ndim
+       call set_coord_limits(i,-halfwidth,halfwidth,'--lim')
+    enddo
+    print "(2(a,1pg10.3),a)",' plot limits: centred box [',-halfwidth,',',halfwidth,']'
+ elseif (have_limits) then
+    !
+    ! set individual column limits from --limits=[min,max,min,max,...]
+    !
+    vals = rflaglist('limits',maxvals,errval=huge(0.),ngot=nvals)
+    if (nvals < 2 .or. mod(nvals,2) /= 0) then
+       print "(a)",' ERROR: --limits requires an even number of comma-separated values'
+       stop
+    endif
+    if (any(vals(1:nvals) >= huge(0.))) then
+       print "(a)",' ERROR: --limits contains a value that is not a number'
+       stop
+    endif
+    if (nvals == 2 .and. ndim >= 1) then
+       !
+       ! if we have exactly 2 values then apply to all coordinate axes (like --lim)
+       !
+       do i=1,ndim
+          call set_coord_limits(i,vals(1),vals(2),'--limits')
+       enddo
+       print "(a,2(1x,1pg10.3))",' plot limits: coordinate axes set to ',vals(1:2)
+    elseif (nvals == 2*ndim .and. ndim >= 1) then
+       ! 
+       ! if we have exactly 2*ndim values then these
+       ! are assumed to apply to coordinate columns
+       !
+       do i=1,ndim
+          call set_coord_limits(i,vals(2*i-1),vals(2*i),'--limits')
+       enddo
+       print "(a,i2,a)",' plot limits: set from --limits for ',ndim,' axes'
+    else
+       !
+       ! otherwise apply to all columns
+       !
+       do icol=1,nvals/2
+          if (icol > numplot) then
+             print "(a,i3)",' ERROR: --limits column index exceeds numplot: ',icol
+             stop
+          endif
+          lim(icol,1:2) = vals(2*icol-1:2*icol)
+       enddo
+       print "(a,i3,a)",' plot limits: overrode ',nvals/2,' columns from --limits'
+    endif
+ else
+    !
+    ! set x,y and z limits from --xmin/--xmax flags
+    !
+    do i=1,ndim
+       if (.not. have_min(i) .and. .not. have_max(i)) cycle
+       call check_coord_column(i,'--'//trim(labelcoord(i,1))//'min/max')
+       xmin = get_command_option(trim(labelcoord(i,1))//'min',default=lim(ix(i),1))
+       xmax = get_command_option(trim(labelcoord(i,1))//'max',default=lim(ix(i),2))
+       lim(ix(i),1:2) = (/xmin,xmax/)
+       print "(a,2(1x,1pg10.3))",' plot limits: '//trim(labelcoord(i,1))//' set to',xmin,xmax
+    enddo
+ endif
+
+end subroutine override_limits_from_flags
+
 end module limits