(lightcurve) added bad pixel fraction to lightcurve outputs, based on…

… number of pixels with tau > 1/3 from individual particles

src/analysis.f90

@@ -421,7 +421,7 @@ subroutine open_analysis(analysistype,required,ncolumns,ndim,ndimV,nsinks)
     endif
     write(headerline,"('#',8(1x,'[',i2.2,1x,a11,']',2x))") &
            1,'time'//trim(labelt),2,'Luminosity',3,'R_{eff}',4,'T_{eff}',&
-           5,'L_{bol}',6,'R_{bb}',7,'T_c'
+           5,'L_{bol}',6,'R_{bb}',7,'T_c',8,'bad pix %'
 
  case('extinction')
     !
@@ -559,7 +559,7 @@ subroutine write_analysis(time,dat,ntot,ntypes,npartoftype,massoftype,&
  real(kind=doub_prec) :: lmin(maxplot),lmax(maxplot),lmean(maxplot),rmsvali
  real(kind=doub_prec), dimension(3) :: xmom,angmom,angmomi,ri,vi
  real                 :: delta,dn,valmin,valmax,valmean,timei
- real                 :: lum,rphoto,tphoto,l_bb,r_bb,t_bb
+ real                 :: lum,rphoto,tphoto,l_bb,r_bb,t_bb,badfrac
  character(len=20)    :: fmtstring
  logical              :: change_coordsys
  real                 :: x0(3),v0(3)
@@ -1260,12 +1260,12 @@ subroutine write_analysis(time,dat,ntot,ntypes,npartoftype,massoftype,&
     return
  case('lightcurve')
     call get_lightcurve(ncolumns,dat,npartoftype,massoftype,iamtype,ndim,ntypes,&
-         lum,rphoto,tphoto,l_bb,r_bb,t_bb,basename(rootname(ifileopen)),ierr)
+         lum,rphoto,tphoto,l_bb,r_bb,t_bb,badfrac,basename(rootname(ifileopen)),ierr)
     print "(4(/,1x,a20,' = ',es9.2))",'Luminosity',lum,'photospheric radius ',rphoto,'photospheric temperature',tphoto
     !
     !--write line to output file
     !
-    if (ierr == 0) write(iunit,"(7(es18.10,1x))") timei,lum,rphoto,tphoto,l_bb,r_bb,t_bb
+    if (ierr == 0) write(iunit,"(8(es18.10,1x))") timei,lum,rphoto,tphoto,l_bb,r_bb,t_bb,badfrac
 
  case('extinction')
     !

src/interpolate3D_opacity.f90

@@ -78,7 +78,7 @@ module interpolate3D_opacity
 
 subroutine interp3D_proj_opacity(x,y,z,pmass,npmass,hh,weight,dat,zorig,itype,npart, &
      xmin,ymin,datsmooth,tausmooth,npixx,npixy,pixwidthx,pixwidthy,zobserver,dscreenfromobserver, &
-     rkappa,zcut,iverbose,exact_rendering,datv,datvpix)
+     rkappa,zcut,iverbose,exact_rendering,datv,datvpix,badpix)
 
  real, parameter :: pi=4.*atan(1.)
  integer, intent(in) :: npart,npixx,npixy,npmass,iverbose
@@ -91,6 +91,8 @@ subroutine interp3D_proj_opacity(x,y,z,pmass,npmass,hh,weight,dat,zorig,itype,np
  ! optional arguments for vector opacity rendering
  real, dimension(:,:),   intent(in),  optional :: datv
  real, dimension(:,:,:), intent(out), optional :: datvpix
+ ! optional argument for checking unresolved pixels in the photosphere
+ real, dimension(npixx,npixy), intent(out), optional :: badpix
 
  integer :: i,ipix,jpix,ipixmin,ipixmax,jpixmin,jpixmax,nused,nsink
  integer, dimension(npart) :: iorder
@@ -116,6 +118,7 @@ subroutine interp3D_proj_opacity(x,y,z,pmass,npmass,hh,weight,dat,zorig,itype,np
  term = 0.
  tausmooth = 0.
  if (present(datvpix)) datvpix = 0.
+ if (present(badpix)) badpix = 0.
  if (pixwidthx <= 0. .or. pixwidthy <= 0) then
     if (iverbose >= -1) print "(1x,a)",'ERROR: pixel width <= 0'
     return
@@ -144,6 +147,7 @@ subroutine interp3D_proj_opacity(x,y,z,pmass,npmass,hh,weight,dat,zorig,itype,np
  !--whether to raytrace backwards from observer, or forwards to observer
  !
  backwards = .false.
+ if (present(badpix)) backwards = .true.
 
  !
  !--setup kernel table if not already set
@@ -370,13 +374,23 @@ subroutine interp3D_proj_opacity(x,y,z,pmass,npmass,hh,weight,dat,zorig,itype,np
                 wab = pixint*dfac
 
                 tau = wab*term
-                fopacity = 1. - exp(-tau)
+                tausmooth(ipix,jpix) = tausmooth(ipix,jpix) + tau
                 !
                 !--render, obscuring previously drawn pixels by relevant amount
                 !  also calculate total optical depth for each pixel
                 !
-                datsmooth(ipix,jpix) = (1.-fopacity)*datsmooth(ipix,jpix) + fopacity*dati
-                tausmooth(ipix,jpix) = tausmooth(ipix,jpix) + tau
+                if (backwards) then
+                   fopacity = exp(-tausmooth(ipix,jpix))*tau
+                   datsmooth(ipix,jpix) = datsmooth(ipix,jpix) + fopacity*dati
+                   if (present(datv)) datvpix(:,ipix,jpix) = datvpix(:,ipix,jpix) + fopacity*datvi(:)
+                   if (present(badpix)) then
+                      if (tausmooth(ipix,jpix) < 1. .and. tau >= 0.33) badpix(ipix,jpix) = 1
+                   endif
+                else
+                   fopacity = 1. - exp(-tau)
+                   datsmooth(ipix,jpix) = (1.-fopacity)*datsmooth(ipix,jpix) + fopacity*dati
+                   if (present(datv)) datvpix(:,ipix,jpix) = (1.-fopacity)*datvpix(:,ipix,jpix) + fopacity*datvi(:)
+                endif
              else
                 !
                 !--SPH kernel - integral through cubic spline
@@ -395,6 +409,9 @@ subroutine interp3D_proj_opacity(x,y,z,pmass,npmass,hh,weight,dat,zorig,itype,np
                       fopacity = exp(-tausmooth(ipix,jpix))*tau
                       datsmooth(ipix,jpix) = datsmooth(ipix,jpix) + fopacity*dati
                       if (present(datv)) datvpix(:,ipix,jpix) = datvpix(:,ipix,jpix) + fopacity*datvi(:)
+                      if (present(badpix)) then
+                         if (tausmooth(ipix,jpix) < 1. .and. tau >= 0.33) badpix(ipix,jpix) = 1
+                      endif
                    else
                       fopacity = 1. - exp(-tau)
                       datsmooth(ipix,jpix) = (1.-fopacity)*datsmooth(ipix,jpix) + fopacity*dati

src/lightcurve.f90

@@ -46,7 +46,7 @@ module lightcurve
 ! Used to generate synthetic lightcurves
 !---------------------------------------------------------
 subroutine get_lightcurve(ncolumns,dat,npartoftype,masstype,itype,ndim,ntypes,&
-                          lum,rphoto,temp,lum_bb,r_bb,Tc,specfile,ierr)
+                          lum,rphoto,temp,lum_bb,r_bb,Tc,badfrac,specfile,ierr)
  use labels,                only:ix,ih,irho,ipmass,itemp,ikappa,ivx,ipmomx
  use limits,                only:lim,get_particle_subset
  use lightcurve_utils,      only:get_temp_from_u
@@ -66,24 +66,24 @@ subroutine get_lightcurve(ncolumns,dat,npartoftype,masstype,itype,ndim,ntypes,&
  use settings_xsecrot,      only:anglex,angley,anglez,taupartdepth
  use rotation,              only:rotate_particles
  use system_utils,          only:get_command_flag
- use readwrite_fits,        only:write_fits_cube,write_fits_image
+ use readwrite_fits,        only:write_fits_cube,write_fits_image,write_fits_image
  integer, intent(in)  :: ncolumns,ntypes,ndim
  integer, intent(in)  :: npartoftype(:)
  integer(kind=int1), intent(in) :: itype(:)
  real,    intent(in)  :: masstype(:)
  real,    intent(in)  :: dat(:,:)
- real,    intent(out) :: lum,rphoto,temp,lum_bb,r_bb,Tc
+ real,    intent(out) :: lum,rphoto,temp,lum_bb,r_bb,Tc,badfrac
  integer, intent(out) :: ierr
  character(len=*), intent(in) :: specfile
  integer :: n,isinktype,npixx,npixy,j,i,k,nfreq
  integer, parameter :: iu1 = 45
  real, dimension(3) :: xmin,xmax
  real, dimension(:),   allocatable :: weight,x,y,z,flux,opacity
  real, dimension(:),   allocatable :: freq,spectrum,bb_spectrum
- real, dimension(:,:), allocatable :: img,taupix,flux_nu,v_on_c
- real, dimension(:,:,:), allocatable :: img_nu
+ 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
- real :: ax,ay,az,xi(3),betaz,lorentz,doppler_factor,doppler_factor_max,taui,lprev
+ real :: ax,ay,az,xi(3),betaz,lorentz,doppler_factor,doppler_factor_max,taui,lprev,badarea
  logical :: relativistic
  character(len=20) :: tmpfile
 
@@ -153,7 +153,7 @@ subroutine get_lightcurve(ncolumns,dat,npartoftype,masstype,itype,ndim,ntypes,&
  !--allocate memory for image
  !
  if (allocated(img) .or. allocated(taupix)) deallocate(img,taupix)
- allocate(img(npixx,npixy),taupix(npixx,npixy))
+ allocate(img(npixx,npixy),taupix(npixx,npixy),badpix(npixx,npixy))
  !
  !--set interpolation weights (w = m/(rho*h^ndim)
  !
@@ -221,10 +221,7 @@ subroutine get_lightcurve(ncolumns,dat,npartoftype,masstype,itype,ndim,ntypes,&
       dat(1:n,ipmass),n,dat(1:n,ih),weight, &
       flux,z,icolourme(1:n), &
       n,xmin(1),xmin(2),img,taupix,npixx,npixy,&
-      dx,dy,zobs,dzobs,opacity,huge(zobs),iverbose,.false.,datv=flux_nu,datvpix=img_nu)
-
- lum = 4.*sum(img)*dx*dy
- print*,'grey luminosity = ',lum,' erg/s'
+      dx,dy,zobs,dzobs,opacity,huge(zobs),iverbose,.false.,datv=flux_nu,datvpix=img_nu,badpix=badpix)
 
  ! integrate flux over all frequencies to give Flux = \int F_\nu d\nu = pi \int B_nu dnu
  do j=1,npixy
@@ -239,7 +236,13 @@ subroutine get_lightcurve(ncolumns,dat,npartoftype,masstype,itype,ndim,ntypes,&
  print "(/,a,2(es10.3,a))",' L_bol = ',lum,' erg/s = ',lum/Lsun,' L_sun'
 
  area = count(taupix >= 1.)*dx*dy
- print "(a,1pg10.3,a)",' emitting area = ',area/au**2,' au^2'
+ badarea = count(badpix > 0.)*dx*dy
+ badfrac = 0.
+ if (area > 0.) badfrac = badarea/area
+ print "(a,1pg10.3,a)",' emitting area = ',area/au**2,' au^2 (pixels where dtau > 1/3)'
+ print "(a,1pg10.3,a,2pf6.2,a)",' unresolved area = ',badarea/au**2,' au^2 (',badfrac,'%)'
+ if (badarea/area > 0.05) print "(/,1x,a,2pf6.2,a)",'WARNING: ',badfrac,'% of photosphere is UNRESOLVED!'
+
  print "(/,a,1pg10.3,a)",' Tmax  = ',(maxval(img)/steboltz)**0.25,' K'
 
  ! effective temperature: total flux equals that of a blackbody at T=Teff
@@ -251,9 +254,12 @@ subroutine get_lightcurve(ncolumns,dat,npartoftype,masstype,itype,ndim,ntypes,&
  ! get integrated spectrum from integrating over all pixels in the image
  if (allocated(spectrum)) deallocate(spectrum,bb_spectrum)
  allocate(spectrum(nfreq),bb_spectrum(nfreq))
+ !$omp parallel do default(none) private(i) &
+ !$omp shared(spectrum,img_nu,npixx,npixy,nfreq,dx,dy)
  do i=1,nfreq
     spectrum(i) = sum(img_nu(i,1:npixx,1:npixy))*dx*dy
  enddo
+ !$omp end parallel do
 
  ! get colour temperature by fitting the blackbody peak
  call get_colour_temperature(spectrum,freq,Tc,freq_max,bb_scale)
@@ -285,9 +291,22 @@ subroutine get_lightcurve(ncolumns,dat,npartoftype,masstype,itype,ndim,ntypes,&
  endif
 
  if (get_command_flag('writefits')) then
-    write(*,*) 'Writing image into FITS cube ...'
-    call write_fits_cube('img_'//trim(specfile)//'.fits',img_nu(:,:,:),(/nfreq,npixx,npixy/),ierr)
+    write(*,"(/,a)") 'Writing total intensity to fits image ...'
+    call write_fits_image('img_'//trim(specfile)//'_mom0.fits',img,(/npixx,npixy/),ierr)
+
+    write(*,"(/,a)") 'Writing image into FITS cube ...'
+    allocate(img_tmp(nfreq,npixx,npixy),stat=ierr)
+    if (ierr == 0) then
+       !$omp parallel do default(none) private(i) &
+       !$omp shared(nfreq,npixx,npixy,img_nu,img_tmp)
+       do i=1,nfreq
+          img_tmp(1:npixx,1:npixy,i) = img_nu(i,1:npixx,1:npixy)
+       enddo
+       !$omp end parallel do
+       call write_fits_cube('img_'//trim(specfile)//'.fits',img_tmp,(/npixx,npixy,nfreq/),ierr)
+    endif
     if (ierr /= 0) write(*,*) 'Error writing to FITS cube !!!'
+    if (allocated(img_tmp)) deallocate(img_tmp)
  endif
 
 end subroutine get_lightcurve