Skip to content
Open
3 changes: 3 additions & 0 deletions docs/lightcurve.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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`.


Expand Down
41 changes: 34 additions & 7 deletions src/lightcurve.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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))
Comment thread
s-neilson marked this conversation as resolved.
Comment thread
coderabbitai[bot] marked this conversation as resolved.

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
Expand All @@ -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))
Expand Down Expand Up @@ -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) &
Comment thread
s-neilson marked this conversation as resolved.
!$omp firstprivate(opacity_factor,doppler_factor) &
!$omp reduction(max:doppler_factor_max)
Expand All @@ -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
Comment thread
s-neilson marked this conversation as resolved.
Comment thread
coderabbitai[bot] marked this conversation as resolved.
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)
Expand Down
Loading