diff --git a/source/merger_trees.operators.export.F90 b/source/merger_trees.operators.export.F90 index 30332409b7..3d884ebf71 100644 --- a/source/merger_trees.operators.export.F90 +++ b/source/merger_trees.operators.export.F90 @@ -263,6 +263,7 @@ subroutine exportOperatePreEvolution(self,tree) if (allocated(snapshotInterpolator)) deallocate(snapshotInterpolator) allocate(snapshotInterpolator) snapshotInterpolator=interpolator(snapshotTime(1:snapshotCount)) + deallocate(snapshotTime) else snapshotCount=0 end if @@ -310,6 +311,7 @@ subroutine exportOperatePreEvolution(self,tree) deallocate(descendantIndex) deallocate(nodeMass ) deallocate(nodeRedshift ) + if (self%snapshotsRequired ) deallocate(nodeSnapshot) if (defaultPositionComponent%positionIsGettable()) deallocate(nodePosition) if (defaultPositionComponent%velocityIsGettable()) deallocate(nodeVelocity) ! Move to the next tree. diff --git a/source/nodes.property_extractor.absolute_magnitude.F90 b/source/nodes.property_extractor.absolute_magnitude.F90 new file mode 100644 index 0000000000..dab4b1eb73 --- /dev/null +++ b/source/nodes.property_extractor.absolute_magnitude.F90 @@ -0,0 +1,217 @@ +!! Copyright 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, +!! 2019, 2020, 2021, 2022, 2023 +!! Andrew Benson +!! +!! This file is part of Galacticus. +!! +!! Galacticus is free software: you can redistribute it and/or modify +!! it under the terms of the GNU General Public License as published by +!! the Free Software Foundation, either version 3 of the License, or +!! (at your option) any later version. +!! +!! Galacticus is distributed in the hope that it will be useful, +!! but WITHOUT ANY WARRANTY; without even the implied warranty of +!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +!! GNU General Public License for more details. +!! +!! You should have received a copy of the GNU General Public License +!! along with Galacticus. If not, see . + +!!{ +Implements a node property extractor class for absolute magnitudes. +!!} + + use :: Galactic_Structure_Options, only : enumerationComponentTypeType + + !![ + + + A node property extractor which extracts stellar absolute magnitudes in all available bands. + + + !!] + type, extends(nodePropertyExtractorTuple) :: nodePropertyExtractorMagnitudesAbsolute + !!{ + A node property extractor which extracts stellar absolute magnitudes in all available bands. + !!} + private + type(enumerationComponentTypeType) :: component + contains + procedure :: elementCount => magnitudesAbsoluteElementCount + procedure :: extract => magnitudesAbsoluteExtract + procedure :: names => magnitudesAbsoluteNames + procedure :: descriptions => magnitudesAbsoluteDescriptions + procedure :: unitsInSI => magnitudesAbsoluteUnitsInSI + end type nodePropertyExtractorMagnitudesAbsolute + + interface nodePropertyExtractorMagnitudesAbsolute + !!{ + Constructors for the ``magnitudesAbsolute'' output analysis class. + !!} + module procedure magnitudesAbsoluteConstructorParameters + module procedure magnitudesAbsoluteConstructorInternal + end interface nodePropertyExtractorMagnitudesAbsolute + +contains + + function magnitudesAbsoluteConstructorParameters(parameters) result(self) + !!{ + Constructor for the {\normalfont \ttfamily magnitudesAbsolute} property extractor class which takes a parameter set as input. + !!} + use :: Input_Parameters , only : inputParameters + use :: Galactic_Structure_Options, only : enumerationComponentTypeEncode + implicit none + type(nodePropertyExtractorMagnitudesAbsolute) :: self + type(inputParameters ), intent(inout) :: parameters + type(varying_string ) :: component + + !![ + + component + parameters + The component from which to extract star formation rate. + + !!] + self=nodePropertyExtractorMagnitudesAbsolute(enumerationComponentTypeEncode(char(component),includesPrefix=.false.)) + !![ + + !!] + return + end function magnitudesAbsoluteConstructorParameters + + function magnitudesAbsoluteConstructorInternal(component) result(self) + !!{ + Internal constructor for the {\normalfont \ttfamily magnitudesAbsolute} property extractor class. + !!} + use :: Input_Parameters, only : inputParameters + implicit none + type(nodePropertyExtractorMagnitudesAbsolute) :: self + type(enumerationComponentTypeType ), intent(in ) :: component + !![ + + !!] + + return + end function magnitudesAbsoluteConstructorInternal + + integer function magnitudesAbsoluteElementCount(self,time) + !!{ + Return the number of elements in the {\normalfont \ttfamily magnitudesAbsolute} property extractor class. + !!} + use :: Stellar_Luminosities_Structure, only : unitStellarLuminosities + implicit none + class (nodePropertyExtractorMagnitudesAbsolute), intent(inout) :: self + double precision , intent(in ) :: time + !$GLC attributes unused :: self + + magnitudesAbsoluteElementCount=unitStellarLuminosities%luminosityOutputCount(time) + return + end function magnitudesAbsoluteElementCount + + function magnitudesAbsoluteExtract(self,node,time,instance) result(magnitudes) + !!{ + Implement a {\normalfont \ttfamily magnitudesAbsolute} property extractor. + !!} + use :: Galacticus_Nodes , only : nodeComponentDisk , nodeComponentSpheroid + use :: Galactic_Structure_Options , only : componentTypeDisk , componentTypeSpheroid + use :: Stellar_Luminosities_Structure, only : stellarLuminosities, unitStellarLuminosities + use :: Error , only : Error_Report + implicit none + double precision , dimension(:) , allocatable :: magnitudes + class (nodePropertyExtractorMagnitudesAbsolute), intent(inout), target :: self + type (treeNode ), intent(inout), target :: node + double precision , intent(in ) :: time + type (multiCounter ), intent(inout), optional :: instance + class (nodeComponentDisk ) , pointer :: disk + class (nodeComponentSpheroid ) , pointer :: spheroid + type (stellarLuminosities ) :: luminosities + integer :: i , j + !$GLC attributes unused :: instance + + allocate(magnitudes(unitStellarLuminosities%luminosityOutputCount(time))) + select case (self%component%ID) + case (componentTypeDisk %ID) + disk => node %disk () + luminosities = disk %luminositiesStellar () + case (componentTypeSpheroid%ID) + spheroid => node %spheroid () + luminosities = spheroid%luminositiesStellar () + case default + luminosities = unitStellarLuminosities + call Error_Report("only 'disk' and 'spheroid' components are supported"//{introspection:location}) + end select + ! Convert luminosities to magnitudes. + j=0 + do i=1,unitStellarLuminosities%luminosityCount() + if (unitStellarLuminosities%isOutput(i,time)) then + j=j+1 + magnitudes(j)=luminosities%luminosity(i) + end if + end do + where (magnitudes > 0.0d0) + ! Where the luminosity is positive, compute the corresponding absolute magnitude. + magnitudes=-2.5d0 & + & *log10(magnitudes) + elsewhere + ! For non-positive luminosities, return the faintest absolute magnitude possible. + magnitudes=+huge (0.0d0 ) + end where + return + end function magnitudesAbsoluteExtract + + subroutine magnitudesAbsoluteNames(self,time,names) + !!{ + Return the names of the {\normalfont \ttfamily magnitudesAbsolute} properties. + !!} + use :: Stellar_Luminosities_Structure, only : unitStellarLuminosities + use :: Galactic_Structure_Options , only : enumerationComponentTypeDecode + implicit none + class (nodePropertyExtractorMagnitudesAbsolute), intent(inout) :: self + double precision , intent(in ) :: time + type (varying_string ), intent(inout), dimension(:) , allocatable :: names + integer :: i , j + !$GLC attributes unused :: self + + allocate(names(unitStellarLuminosities%luminosityOutputCount(time))) + j=0 + do i=1,unitStellarLuminosities%luminosityCount() + if (unitStellarLuminosities%isOutput(i,time)) then + j=j+1 + names(j)=enumerationComponentTypeDecode(self%component)//'MagnitudeAbsoluteStellar:'//unitStellarLuminosities%name(i) + end if + end do + return + end subroutine magnitudesAbsoluteNames + + subroutine magnitudesAbsoluteDescriptions(self,time,descriptions) + !!{ + Return descriptions of the {\normalfont \ttfamily magnitudesAbsolute} property extractor class. + !!} + use :: Stellar_Luminosities_Structure, only : unitStellarLuminosities + implicit none + class (nodePropertyExtractorMagnitudesAbsolute), intent(inout) :: self + double precision , intent(in ) :: time + type (varying_string ), intent(inout), dimension(:) , allocatable :: descriptions + !$GLC attributes unused :: self + + allocate(descriptions(unitStellarLuminosities%luminosityOutputCount(time))) + descriptions=var_str('Stellar absolute magnitude [AB system].') + return + end subroutine magnitudesAbsoluteDescriptions + + function magnitudesAbsoluteUnitsInSI(self,time) result(unitsInSI) + !!{ + Return the units of absolute magnitude in the SI system. + !!} + use :: Stellar_Luminosities_Structure, only : unitStellarLuminosities + implicit none + double precision , allocatable , dimension(:) :: unitsInSI + class (nodePropertyExtractorMagnitudesAbsolute), intent(inout) :: self + double precision , intent(in ) :: time + !$GLC attributes unused :: self + + allocate(unitsInSI(unitStellarLuminosities%luminosityOutputCount(time))) + unitsInSI=0.0d0 + return + end function magnitudesAbsoluteUnitsInSI + diff --git a/source/nodes.property_extractor.apparent_magnitude.F90 b/source/nodes.property_extractor.apparent_magnitude.F90 new file mode 100644 index 0000000000..ae72a3ea9d --- /dev/null +++ b/source/nodes.property_extractor.apparent_magnitude.F90 @@ -0,0 +1,249 @@ +!! Copyright 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, +!! 2019, 2020, 2021, 2022, 2023 +!! Andrew Benson +!! +!! This file is part of Galacticus. +!! +!! Galacticus is free software: you can redistribute it and/or modify +!! it under the terms of the GNU General Public License as published by +!! the Free Software Foundation, either version 3 of the License, or +!! (at your option) any later version. +!! +!! Galacticus is distributed in the hope that it will be useful, +!! but WITHOUT ANY WARRANTY; without even the implied warranty of +!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +!! GNU General Public License for more details. +!! +!! You should have received a copy of the GNU General Public License +!! along with Galacticus. If not, see . + +!!{ +Implements a node property extractor class for apparent magnitudes. +!!} + + use :: Galactic_Structure_Options, only : enumerationComponentTypeType + use :: Cosmology_Functions , only : cosmologyFunctionsClass + + !![ + + + A node property extractor which extracts stellar apparent magnitudes in all available bands. + + + !!] + type, extends(nodePropertyExtractorTuple) :: nodePropertyExtractorMagnitudesApparent + !!{ + A node property extractor which extracts stellar apparent magnitudes in all available bands. + !!} + private + class(cosmologyFunctionsClass ), pointer :: cosmologyFunctions_ => null() + type (enumerationComponentTypeType) :: component + contains + final :: magnitudesApparentDestructor + procedure :: elementCount => magnitudesApparentElementCount + procedure :: extract => magnitudesApparentExtract + procedure :: names => magnitudesApparentNames + procedure :: descriptions => magnitudesApparentDescriptions + procedure :: unitsInSI => magnitudesApparentUnitsInSI + end type nodePropertyExtractorMagnitudesApparent + + interface nodePropertyExtractorMagnitudesApparent + !!{ + Constructors for the ``magnitudesApparent'' output analysis class. + !!} + module procedure magnitudesApparentConstructorParameters + module procedure magnitudesApparentConstructorInternal + end interface nodePropertyExtractorMagnitudesApparent + +contains + + function magnitudesApparentConstructorParameters(parameters) result(self) + !!{ + Constructor for the {\normalfont \ttfamily magnitudesApparent} property extractor class which takes a parameter set as input. + !!} + use :: Input_Parameters , only : inputParameters + use :: Galactic_Structure_Options, only : enumerationComponentTypeEncode + implicit none + type (nodePropertyExtractorMagnitudesApparent) :: self + type (inputParameters ), intent(inout) :: parameters + class(cosmologyFunctionsClass ), pointer :: cosmologyFunctions_ + type (varying_string ) :: component + + !![ + + component + parameters + The component from which to extract star formation rate. + + + !!] + self=nodePropertyExtractorMagnitudesApparent(enumerationComponentTypeEncode(char(component),includesPrefix=.false.),cosmologyFunctions_) + !![ + + + !!] + return + end function magnitudesApparentConstructorParameters + + function magnitudesApparentConstructorInternal(component,cosmologyFunctions_) result(self) + !!{ + Internal constructor for the {\normalfont \ttfamily magnitudesApparent} property extractor class. + !!} + use :: Input_Parameters, only : inputParameters + implicit none + type (nodePropertyExtractorMagnitudesApparent) :: self + type (enumerationComponentTypeType ), intent(in ) :: component + class(cosmologyFunctionsClass ), intent(in ), target :: cosmologyFunctions_ + !![ + + !!] + + return + end function magnitudesApparentConstructorInternal + + subroutine magnitudesApparentDestructor(self) + !!{ + Destructor for the {\normalfont \ttfamily magnitudesApparent} property extractor class. + !!} + implicit none + type(nodePropertyExtractorMagnitudesApparent), intent(inout) :: self + + !![ + + !!] + return + end subroutine magnitudesApparentDestructor + + integer function magnitudesApparentElementCount(self,time) + !!{ + Return the number of elements in the {\normalfont \ttfamily magnitudesApparent} property extractor class. + !!} + use :: Stellar_Luminosities_Structure, only : unitStellarLuminosities + implicit none + class (nodePropertyExtractorMagnitudesApparent), intent(inout) :: self + double precision , intent(in ) :: time + !$GLC attributes unused :: self + + magnitudesApparentElementCount=unitStellarLuminosities%luminosityOutputCount(time) + return + end function magnitudesApparentElementCount + + function magnitudesApparentExtract(self,node,time,instance) result(magnitudes) + !!{ + Implement a {\normalfont \ttfamily magnitudesApparent} property extractor. + !!} + use :: Galacticus_Nodes , only : nodeComponentDisk , nodeComponentSpheroid + use :: Galactic_Structure_Options , only : componentTypeDisk , componentTypeSpheroid + use :: Stellar_Luminosities_Structure, only : stellarLuminosities, unitStellarLuminosities + use :: Error , only : Error_Report + implicit none + double precision , dimension(:) , allocatable :: magnitudes + class (nodePropertyExtractorMagnitudesApparent), intent(inout), target :: self + type (treeNode ), intent(inout), target :: node + double precision , intent(in ) :: time + type (multiCounter ), intent(inout), optional :: instance + class (nodeComponentDisk ) , pointer :: disk + class (nodeComponentSpheroid ) , pointer :: spheroid + type (stellarLuminosities ) :: luminosities + integer :: i , j + double precision :: distanceModulus, distanceLuminosity + !$GLC attributes unused :: instance + + allocate(magnitudes(unitStellarLuminosities%luminosityOutputCount(time))) + select case (self%component%ID) + case (componentTypeDisk %ID) + disk => node %disk () + luminosities = disk %luminositiesStellar () + case (componentTypeSpheroid%ID) + spheroid => node %spheroid () + luminosities = spheroid%luminositiesStellar () + case default + luminosities = unitStellarLuminosities + call Error_Report("only 'disk' and 'spheroid' components are supported"//{introspection:location}) + end select + ! Find the distance modulus, including the extra -2.5log₁₀(1+z) term that accounts for compression of photon frequencies due to redshifting. + distanceLuminosity=self%cosmologyFunctions_%distanceLuminosity(time) + if (distanceLuminosity > 0.0d0) then + distanceModulus=+25.0d0 & + & + 5.0d0*log10( distanceLuminosity ) & + & + 2.5d0*log10(self%cosmologyFunctions_%expansionFactor (time)) + ! Convert luminosities to magnitudes. + j=0 + do i=1,unitStellarLuminosities%luminosityCount() + if (unitStellarLuminosities%isOutput(i,time)) then + j=j+1 + magnitudes(j)=luminosities%luminosity(i) + end if + end do + where (magnitudes > 0.0d0) + ! Where the luminosity is positive, compute the corresponding apparent magnitude. + magnitudes=-2.5d0 & + & *log10(magnitudes) & + & +distanceModulus + elsewhere + ! For non-positive luminosities, return the faintest apparent magnitude possible. + magnitudes=+huge (0.0d0 ) + end where + else + ! Luminosity distance is zero, so set apparent magnitudes to the brightest magnitude possible. + magnitudes =-huge (0.0d0 ) + end if + return + end function magnitudesApparentExtract + + subroutine magnitudesApparentNames(self,time,names) + !!{ + Return the names of the {\normalfont \ttfamily magnitudesApparent} properties. + !!} + use :: Stellar_Luminosities_Structure, only : unitStellarLuminosities + use :: Galactic_Structure_Options , only : enumerationComponentTypeDecode + implicit none + class (nodePropertyExtractorMagnitudesApparent), intent(inout) :: self + double precision , intent(in ) :: time + type (varying_string ), intent(inout), dimension(:) , allocatable :: names + integer :: i , j + !$GLC attributes unused :: self + + allocate(names(unitStellarLuminosities%luminosityOutputCount(time))) + j=0 + do i=1,unitStellarLuminosities%luminosityCount() + if (unitStellarLuminosities%isOutput(i,time)) then + j=j+1 + names(j)=enumerationComponentTypeDecode(self%component)//'MagnitudeApparentStellar:'//unitStellarLuminosities%name(i) + end if + end do + return + end subroutine magnitudesApparentNames + + subroutine magnitudesApparentDescriptions(self,time,descriptions) + !!{ + Return descriptions of the {\normalfont \ttfamily magnitudesApparent} property extractor class. + !!} + use :: Stellar_Luminosities_Structure, only : unitStellarLuminosities + implicit none + class (nodePropertyExtractorMagnitudesApparent), intent(inout) :: self + double precision , intent(in ) :: time + type (varying_string ), intent(inout), dimension(:) , allocatable :: descriptions + !$GLC attributes unused :: self + + allocate(descriptions(unitStellarLuminosities%luminosityOutputCount(time))) + descriptions=var_str('Stellar apparent magnitude [AB system].') + return + end subroutine magnitudesApparentDescriptions + + function magnitudesApparentUnitsInSI(self,time) result(unitsInSI) + !!{ + Return the units of apparent magnitude in the SI system. + !!} + use :: Stellar_Luminosities_Structure, only : unitStellarLuminosities + implicit none + double precision , allocatable , dimension(:) :: unitsInSI + class (nodePropertyExtractorMagnitudesApparent), intent(inout) :: self + double precision , intent(in ) :: time + !$GLC attributes unused :: self + + allocate(unitsInSI(unitStellarLuminosities%luminosityOutputCount(time))) + unitsInSI=0.0d0 + return + end function magnitudesApparentUnitsInSI + diff --git a/source/objects.merger_tree_data.F90 b/source/objects.merger_tree_data.F90 index be7ba61788..eed51f3653 100644 --- a/source/objects.merger_tree_data.F90 +++ b/source/objects.merger_tree_data.F90 @@ -1588,7 +1588,7 @@ subroutine Merger_Tree_Data_Structure_Export_Galacticus(mergerTrees,outputFileNa if (mergerTrees%hasNodeMass ) call forestHalos%writeDataset(mergerTrees%nodeMass ,"nodeMass" ,"The mass of each node." ,appendTo=appendActual ) if (mergerTrees%hasRedshift ) call forestHalos%writeDataset(mergerTrees%redshift ,"redshift" ,"The redshift of each node." ,appendTo=appendActual ) if (mergerTrees%hasNodeMass200Mean ) call forestHalos%writeDataset(mergerTrees%nodeMass200Mean ,"nodeMass200Mean" ,"The M200 mass of each node (200 * mean density)." ,appendTo=appendActual ) - if (mergerTrees%hasNodeMass200Crit ) call forestHalos%writeDataset(mergerTrees%nodeMass200Crit ,"nodeMass200Crit" ,"The M200 mass of each node (200 * crit density(." ,appendTo=appendActual ) + if (mergerTrees%hasNodeMass200Crit ) call forestHalos%writeDataset(mergerTrees%nodeMass200Crit ,"nodeMass200Crit" ,"The M200 mass of each node (200 * crit density)." ,appendTo=appendActual ) if (mergerTrees%hasScaleFactor ) call forestHalos%writeDataset(mergerTrees%scaleFactor ,"scaleFactor" ,"The scale factor of each node." ,appendTo=appendActual ) if (mergerTrees%hasPositionX ) call forestHalos%writeDataset(mergerTrees%position ,"position" ,"The position of each node." ,appendTo=appendActual,appendDimension=2) if (mergerTrees%hasVelocityX ) call forestHalos%writeDataset(mergerTrees%velocity ,"velocity" ,"The velocity of each node." ,appendTo=appendActual,appendDimension=2) @@ -1974,8 +1974,8 @@ subroutine Merger_Tree_Data_Structure_Export_IRATE(mergerTrees,outputFileName,hd ! Output merger tree datasets. call mergerTreesGroup%writeDataset(mergerTrees%snapshot ,"HaloSnapshot" ,"The snapshot of each halo." ,appendTo=appendActual) call mergerTreesGroup%writeDataset(mergerTrees%nodeIndex ,"HaloID" ,"The index of each halo." ,appendTo=appendActual) - call mergerTreesGroup%writeDataset(mergerTrees%descendantIndex ,"DescendantID" ,"The index of each descendant halo." ,appendTo=appendActual) - call mergerTreesGroup%writeDataset( descendantSnapshot,"DescendantSnapshot","The snapshot of each descendant halo.",appendTo=appendActual) + call mergerTreesGroup%writeDataset(mergerTrees%descendantIndex ,"DescendentID" ,"The index of each descendant halo." ,appendTo=appendActual) + call mergerTreesGroup%writeDataset( descendantSnapshot,"DescendentSnapshot","The snapshot of each descendant halo.",appendTo=appendActual) if (mergerTrees%hasHostIndex) call mergerTreesGroup%writeDataset(mergerTrees%hostIndex ,"HostID" ,"The index of each host halo." ,appendTo=appendActual) call mergerTreesGroup%writeDataset(mergerTrees%treeNodeCount ,"HalosPerTree" ,"Number of halos in each tree." ,appendTo=appendActual) call mergerTreesGroup%writeDataset(mergerTrees%forestID ,"TreeID" ,"Unique index of tree." ,appendTo=appendActual)