diff --git a/docs/lightcurve.rst b/docs/lightcurve.rst index 6f345173..805f5550 100755 --- a/docs/lightcurve.rst +++ b/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 is by default equal to the extent of the particle positions in the dump file. Any number of these automatically calculated bounds can be overriden with the ``SPLASH_MARGIN_XMIN``, ``SPLASH_MARGIN_XMAX``, ``SPLASH_MARGIN_YMIN`` and ``SPLASH_MARGIN_YMAX`` environment variables or the ``--xmin``, ``--xmax``, ``--ymin`` and ``--ymax`` flags. The units of these bounds are the same as the units of the particle positions. 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 (called the spectral hardening factor) and constant offset. The scaling factor is set using the flag ``--fcol`` and the constant offset with the flag ``--ocol``. The scaling factor is applied first and the constant offset is added to the scaled temperature. + Opacity :math:`\kappa` used in the ray trace is described in :ref:`sec:lckappa`. diff --git a/src/lightcurve.f90 b/src/lightcurve.f90 index 472d1251..66a8f4a2 100644 --- a/src/lightcurve.f90 +++ b/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,27 @@ 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') + print "(a,f8.2)"," Temperature offset: ",o_col + + !If the user does not provide a limit for the render region it is calculated automatically. xmin(1:ndim) = lim(ix(1:ndim),1) - xmax(1:ndim) = lim(ix(1:ndim),2) + xmax(1:ndim) = lim(ix(1:ndim),2) + xmin(1) = renvironment('SPLASH_MARGIN_XMIN', xmin(1)) + xmax(1) = renvironment('SPLASH_MARGIN_XMAX', xmax(1)) + xmin(2) = renvironment('SPLASH_MARGIN_YMIN', xmin(2)) + xmax(2) = renvironment('SPLASH_MARGIN_YMAX', xmax(2)) + + if (xmax(1) <= xmin(1)) then + print *,"ERROR: maximum x axis render bound less than or equal to minimum x axis render region bound, aborting..." + ierr = 3 + return + endif + if (xmax(2) <= xmin(2)) then + print *,"ERROR: maximum y axis render bound less than or equal to minimum y axis render region bound, aborting..." + ierr = 4 + return + endif ! !--set number of particles to use in the interpolation routines ! and allocate memory for weights @@ -137,7 +156,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 = 5 return endif x(1:n) = dat(1:n,ix(1)) @@ -213,8 +232,8 @@ subroutine get_lightcurve(ncolumns,dat,npartoftype,masstype,itype,ndim,ntypes,& opacity_factor = 1. 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(n,nfreq,freq,flux,flux_nu,dat,h,weight,radkernel,ierr) & + !$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 +258,20 @@ 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 + if (tempi <= 0.0) ierr = 6 + 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 (ierr /= 6) flux_nu(:,i) = B_nu(tempi,freq*doppler_factor) enddo !$omp end parallel do + + if (ierr == 6) then + print *,"ERROR: choice of temperature offset has resulted in a zero or negative temperature for at least one particle, aborting..." + return + endif + if (relativistic) print*,' max relativistic correction=',doppler_factor_max if (allocated(img_nu)) deallocate(img_nu)