MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PSD.cc
Go to the documentation of this file.
1 //Copyright (c) 2013-2020, The MercuryDPM Developers Team. All rights reserved.
2 //For the list of developers, see <http://www.MercuryDPM.org/Team>.
3 //
4 //Redistribution and use in source and binary forms, with or without
5 //modification, are permitted provided that the following conditions are met:
6 // * Redistributions of source code must retain the above copyright
7 // notice, this list of conditions and the following disclaimer.
8 // * Redistributions in binary form must reproduce the above copyright
9 // notice, this list of conditions and the following disclaimer in the
10 // documentation and/or other materials provided with the distribution.
11 // * Neither the name MercuryDPM nor the
12 // names of its contributors may be used to endorse or promote products
13 // derived from this software without specific prior written permission.
14 //
15 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16 //ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17 //WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18 //DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
19 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21 //LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22 //ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23 //(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24 //SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 
26 #include <DPMBase.h>
27 #include "PSD.h"
28 
33 {
34  for (auto& momenta: momenta_)
35  momenta = 0;
36 }
37 
41 PSD::PSD(const PSD& other)
42 {
44  momenta_ = other.momenta_;
45 }
46 
51 PSD::~PSD()
52 = default;
53 
54 
59 PSD* PSD::copy() const
60 {
61 #ifdef DEBUG_CONSTRUCTOR
62  std::cout << "PSD::copy() const finished" << std::endl;
63 #endif
64  return new PSD(*this);
65 }
66 
71 void PSD::printPSD() const
72 {
73  if (particleSizeDistribution_.front().radius > 1e-1)
74  {
75  for (const auto p : particleSizeDistribution_)
76  {
77  logger(INFO, "%m\t %\%", p.radius, p.probability * 100);
78  }
79  }
80  else if (particleSizeDistribution_.front().radius > 1e-4)
81  {
82  for (const auto p : particleSizeDistribution_)
83  {
84  logger(INFO, "%mm\t %\%", p.radius * 1000, p.probability * 100);
85  }
86  }
87  else if (particleSizeDistribution_.front().radius > 1e-7)
88  {
89  for (const auto p : particleSizeDistribution_)
90  {
91  logger(INFO, "%um\t %\%", p.radius * 1e6, p.probability * 100);
92  }
93  }
94  else
95  {
96  for (const auto p : particleSizeDistribution_)
97  {
98  logger(INFO, "%nm\t %\%", p.radius * 1e9, p.probability * 100);
99  }
100  }
101 }
102 
110 {
111  // draw a number between 0 and 1, uniformly distributed
113  static std::mt19937 gen(0);
114  static std::uniform_real_distribution<Mdouble> dist(0, 1);
115  Mdouble prob = dist(gen);
116  // find the interval [low,high] of the psd in which the number sits
117  auto high = std::lower_bound(particleSizeDistribution_.begin(), particleSizeDistribution_.end(), prob);
118  auto low = std::max(particleSizeDistribution_.begin(), high - 1);
119  Mdouble rMin = low->radius;
120  Mdouble rMax = high->radius;
121  Mdouble pMin = low->probability;
122  Mdouble pMax = high->probability;
123  // interpolate linearly between [low.radius,high.radius] (assumption: CDF is piecewise linear)
124  Mdouble a = (prob - pMin) / (pMax - pMin);
125  return a * rMin + (1 - a) * rMax;
126 }
127 
138 {
139  // initialize the particleNumberPerClass vector if empty.
140  if (nParticlesPerClass_.empty())
142  // initialize the particleNumberPerClass vector if empty.
143  if (volumePerClass_.empty())
145 
146  for (auto it = particleSizeDistribution_.end() - 1; it != particleSizeDistribution_.begin(); --it)
147  {
148  static std::mt19937 gen(0);
149  static std::uniform_real_distribution<Mdouble> dist(0, 1);
150  Mdouble prob = dist(gen);
151  // map the probability to the current class interval
152  prob = (it - 1)->probability + (it->probability - (it - 1)->probability) * prob;
153  // find the interval [low,high] of the psd in which the number sits
154  auto high = std::lower_bound(particleSizeDistribution_.begin(), particleSizeDistribution_.end(), prob);
155  auto low = std::max(particleSizeDistribution_.begin(), high - 1);
156  Mdouble rMin = low->radius;
157  Mdouble rMax = high->radius;
158  Mdouble pMin = low->probability;
159  Mdouble pMax = high->probability;
160  // interpolate linearly between [low.radius,high.radius] (assumption: CDF is piecewise linear)
161  int index = std::distance(particleSizeDistribution_.begin(), high);
162  Mdouble a = (prob - pMin) / (pMax - pMin);
163  Mdouble rad = a * rMin + (1 - a) * rMax;
164  // Compare inserted volume against volumeAllowed
165  Mdouble radVolume = 4.0 / 3.0 * constants::pi * rad * rad * rad;
168  Mdouble volumeAllowed = particleSizeDistribution_[index].probability * volume;
171  if (volumePerClass_[index] > volumeAllowed)
172  {
173  continue;
174  }
175  else if (radVolume + volumePerClass_[index] > volumeAllowed)
176  {
177  Mdouble differenceVolumeLow = -(volumePerClass_[index] - volumeAllowed);
178  Mdouble differenceVolumeHigh = radVolume + volumePerClass_[index] - volumeAllowed;
179  Mdouble volumeRatio = differenceVolumeLow / differenceVolumeHigh;
180  // If volumeRatio > 1 it will insert anyways, because it should be no problem for the distribution
181  prob = dist(gen);
182  if (prob <= volumeRatio)
183  {
184  ++nParticlesPerClass_[index];
185  volumePerClass_[index] += radVolume;
186  return rad;
187  }
188  else
189  {
190  continue;
191  }
192  }
193  else
194  {
195  ++nParticlesPerClass_[index];
196  volumePerClass_[index] += radVolume;
197  return rad;
198  }
199  }
200  return 0;
201 }
202 
203 
211 {
212  // ensure interval of probabilities to be [0,1]
213  for (auto p = 0; p < particleSizeDistribution_.size(); ++p)
214  particleSizeDistribution_[p].probability =
215  particleSizeDistribution_[p].probability / particleSizeDistribution_.back().probability;
216  // check whether the distribution is cumulative
217  for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
218  {
219  logger.assert_always(it->probability >= (it - 1)->probability, "psd is not cumulative", true);
220  }
221  // cdf needs to start with a probability of zero
222  if (particleSizeDistribution_[0].probability != 0)
223  {
224  logger(INFO, "adding a zero at the beginning of the psd");
226  {std::max(1e-3 * particleSizeDistribution_[0].radius,
227  2 * particleSizeDistribution_[0].radius -
228  particleSizeDistribution_[1].radius), 0});
229  }
230  // cdf needs to end with a probability of one
231  if (particleSizeDistribution_.back().probability < 1)
232  {
233  logger(INFO, "adding a one at the end of the psd");
234  particleSizeDistribution_.push_back({particleSizeDistribution_.back().radius, 1});
235  }
236  // cut off equal subsequent values on the end of the cdf to remove zero size classes.
237  for (auto i = particleSizeDistribution_.end() - 1; i != particleSizeDistribution_.begin(); i--)
238  {
239  if (i->probability == (i - 1)->probability)
240  {
242  .end(), i));
243  }
244  }
245 }
246 
252 {
253  Mdouble sum = 0;
254  // check if the psd starts with zero probability (i.e. p_0=0)
255  if (particleSizeDistribution_[0].probability != 0)
256  {
257  logger(INFO, "adding a zero at the beginning of PDF");
259  {std::max(1e-3 * particleSizeDistribution_[0].radius,
260  2 * particleSizeDistribution_[0].radius -
261  particleSizeDistribution_[1].radius), 0});
262  }
263  // sum up probabilities to check if it equals to unity (sum(p_i)=1)
264  for (auto& it : particleSizeDistribution_)
265  {
266  sum += it.probability;
267  }
268  // ensure interval of probabilities to be [0,1] by normalization.
269  for (auto& it : particleSizeDistribution_)
270  it.probability /= sum;
271  // if the sum of probabilities is not equal to unity, sum up the normalized probabilities again
272  if (sum - 1 > 1e-6)
273  {
274  sum = 0;
275  for (auto& it: particleSizeDistribution_)
276  {
277  sum += it.probability;
278  }
279  }
280  logger.assert_debug(sum - 1 < 1e-6, "PDF is not valid: Integral of PDF is not equal to unity");
281 }
282 
290 void PSD::setDistributionUniform(Mdouble radMin, Mdouble radMax, int numberOfBins)
291 {
292  if (!particleSizeDistribution_.empty())
293  {
295  }
296  std::vector<Mdouble> probabilities = helpers::linspace(0.0, 1.0, numberOfBins);
297  std::vector<Mdouble> radii = helpers::linspace(radMin, radMax, numberOfBins);
298  // combine radii and probabilities
299  for (int i = 0; i < radii.size(); ++i)
300  {
301  particleSizeDistribution_.push_back({radii[i], probabilities[i]});
302  }
304 }
305 
314 void PSD::setDistributionNormal(Mdouble mean, Mdouble standardDeviation, int numberOfBins)
315 {
316  logger.assert_always(mean > 3 * standardDeviation,
317  "Reduce standardDeviation of your normal distribution to avoid negative radii; The mean "
318  "should be greater than 3*standardDeviation");
319  if (!particleSizeDistribution_.empty())
320  {
322  }
323  Mdouble radMin = mean - 3 * standardDeviation;
324  Mdouble radMax = mean + 3 * standardDeviation;
325  std::vector<Mdouble> radii = helpers::linspace(radMin, radMax, numberOfBins);
326  std::vector<Mdouble> probabilities;
327  for (int i = 0; i < radii.size(); i++)
328  {
329  Mdouble probability = 0.5 * (1 + erf((radii[i] - mean) / (sqrt(2) * standardDeviation)));
330  probabilities.push_back(probability);
331  }
332  for (int j = 0; j < radii.size(); j++)
333  {
334  particleSizeDistribution_.push_back({radii[j], probabilities[j]});
335  }
337 }
338 
352 void PSD::setPSDFromVector(std::vector<RadiusAndProbability> psdVector, TYPE PSDType)
353 {
354  particleSizeDistribution_ = psdVector;
355  switch (PSDType)
356  {
357  // default case is a cumulative number distribution (CND)
358  default:
360  break;
367  break;
374  break;
380  break;
385  break;
389  break;
394  break;
399  break;
400  }
401 }
402 
418 void PSD::setPSDFromCSV(const std::string& fileName, TYPE PSDType, bool headings, Mdouble unitScalingFactorRadii)
419 {
420  // read in csv file using the csvReader class
421  csvReader csv;
422  csv.setHeader(headings);
423  csv.read(fileName);
424  std::vector<Mdouble> radii = csv.getFirstColumn(unitScalingFactorRadii);
425  std::vector<Mdouble> probabilities = csv.getSecondColumn(1.0);
426  logger.assert_always(radii.size() == probabilities.size(), "The radii and probabilities vector have to be the "
427  "same size");
428  // combine radii and probabilities
429  for (int i = 0; i < radii.size(); ++i)
430  {
431  particleSizeDistribution_.push_back({radii[i], probabilities[i]});
432  }
433  logger.assert_always(!particleSizeDistribution_.empty(), "PSD cannot be empty");
434  switch (PSDType)
435  {
436  // default case is a CUMULATIVE_NUMBER_DISTRIBUTION
437  default:
439  break;
446  break;
452  break;
459  break;
463  break;
468  break;
473  break;
478  break;
479  }
480  logger(INFO, "A PSD with size ratio % is now ready to be set", getSizeRatio());
481 }
482 
488 {
489  // add up probabilities
490  for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
491  {
492  it->probability = std::min(1.0, it->probability + (it - 1)->probability);
493  }
494  // check whether cumulative distribution is valid
496 }
497 
503 {
504  // subtract cumulative probabilities
505  Mdouble probabilityOld = 0, probabilityCDF;
506  for (auto& it : particleSizeDistribution_)
507  {
508  probabilityCDF = it.probability;
509  it.probability -= probabilityOld;
510  probabilityOld = probabilityCDF;
511  }
512  // check whether probability density distribution is valid
514 }
515 
523 {
524  Mdouble sum = 0;
525  switch (PDFType)
526  {
527  default:
528  logger(ERROR, "Wrong PDFType");
529  break;
531  for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
532  {
533  it->probability /= 0.5 * (it->radius + (it - 1)->radius);
534  // old conversion
535  //p.probability /= p.radius;
536  // sum up probabilities
537  sum += it->probability;
538  }
539  // normalize
540  for (auto& p : particleSizeDistribution_)
541  {
542  p.probability /= sum;
543  }
544  break;
546  for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
547  {
548  it->probability /=
549  1.0 / 3.0 * (square(it->radius) + it->radius * (it - 1)->radius + square((it - 1)->radius));
550  // old conversion
551  //p.probability /= square(p.radius);
552  // sum up probabilities
553  sum += it->probability;
554  }
555  // normalize
556  for (auto& p : particleSizeDistribution_)
557  {
558  p.probability /= sum;
559  }
560  break;
562  for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
563  {
564  it->probability /=
565  0.25 * (square(it->radius) + square((it - 1)->radius)) * (it->radius + (it - 1)->radius);
566  // old conversion
567  //p.probability /= cubic(p.radius);
568  // sum up probabilities
569  sum += it->probability;
570  }
571  // normalize
572  for (auto& p : particleSizeDistribution_)
573  {
574  p.probability /= sum;
575  }
576  break;
577  }
578 }
579 
585 {
586  Mdouble sum = 0;
587  for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
588  {
589  it->probability *=
590  0.25 * (square(it->radius) + square((it - 1)->radius)) * (it->radius + (it - 1)->radius);
591  // old conversion
592  //p.probability /= cubic(p.radius);
593  // sum up probabilities
594  sum += it->probability;
595  }
596  // normalize
597  for (auto& p : particleSizeDistribution_)
598  {
599  p.probability /= sum;
600  }
601 }
602 
610 {
611  Mdouble sum = 0;
612  switch (CDFType)
613  {
614  default:
615  logger(ERROR, "Wrong CDFType");
616  break;
618  for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
619  {
620  // add conversion here
621 
622  // sum up probabilities
623  sum += it->probability;
624  }
625  // normalize
626  for (auto& p : particleSizeDistribution_)
627  {
628  p.probability /= sum;
629  }
630  break;
632  for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
633  {
634  // add conversion here
635 
636  // sum up probabilities
637  sum += it->probability;
638  }
639  // normalize
640  for (auto& p : particleSizeDistribution_)
641  {
642  p.probability /= sum;
643  }
644  break;
646  for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
647  {
648  // add conversion here
649 
650  // sum up probabilities
651  sum += it->probability;
652  }
653  // normalize
654  for (auto& p : particleSizeDistribution_)
655  {
656  p.probability /= sum;
657  }
658  break;
659  }
660 }
661 
668 void PSD::cutoffCumulativeNumber(Mdouble quantileMin, Mdouble quantileMax, Mdouble minPolydispersity)
669 {
670  Mdouble radiusMin = getRadiusByQuantile(quantileMin);
671  Mdouble radiusMax = getRadiusByQuantile(quantileMax);
672  // to get a minimum polydispersity at the base
673  Mdouble radiusMinCut = std::min(radiusMin * (1 + minPolydispersity), radiusMax);
674  // cut off min
675  while (particleSizeDistribution_.front().radius <= radiusMinCut)
676  {
677  particleSizeDistribution_.erase(particleSizeDistribution_.begin());
678  }
679  particleSizeDistribution_.insert(particleSizeDistribution_.begin(), {radiusMinCut, quantileMin});
680  particleSizeDistribution_.insert(particleSizeDistribution_.begin(), {radiusMin, 0});
681  // cut off max
682  while (particleSizeDistribution_.back().radius >= radiusMax)
683  {
684  particleSizeDistribution_.pop_back();
685  }
686  Mdouble radiusMaxCut = std::max(radiusMax - (1 - minPolydispersity) * (radiusMax - particleSizeDistribution_.back()
687  .radius), radiusMin);
688  particleSizeDistribution_.push_back({radiusMaxCut, quantileMax});
689  particleSizeDistribution_.push_back({radiusMax, 1});
690  logger(INFO, "PSD was cut successfully and now has a size ratio of %", getSizeRatio());
691 }
692 
702 void PSD::cutoffAndSqueezeCumulative(Mdouble quantileMin, Mdouble quantileMax, Mdouble squeeze, Mdouble
703 minPolydispersity)
704 {
705  Mdouble r50 = 0.5 * PSD::getNumberDx(50);
706  // cut off
707  cutoffCumulativeNumber(quantileMin, quantileMax, minPolydispersity);
708  // squeeze psd
709  for (auto& p: particleSizeDistribution_)
710  {
711  p.radius = r50 + (p.radius - r50) * squeeze;
712  }
713 }
714 
721 {
722  return 2.0 * getRadiusByQuantile(x / 100);
723 }
724 
731 {
732  PSD psd = *this;
736  return 2.0 * psd.getRadiusByQuantile(x / 100);
737 }
738 
745 {
746  logger.assert_always(quantile <= 1 && quantile >= 0, "quantile is not between 0 and 1");
747  // find the quantile corresponding to the PSD
748  auto high = std::lower_bound(particleSizeDistribution_.begin(), particleSizeDistribution_.end(), quantile);
749  auto low = std::max(particleSizeDistribution_.begin(), high - 1);
750  if (high->probability == low->probability)
751  return high->radius;
752  else
753  return low->radius + (high->radius - low->radius) * (quantile - low->probability) /
754  (high->probability - low->probability);
755 }
756 
763 {
764  PSD psd = *this;
767  Mdouble mean = 0;
768  for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
769  mean += it->probability * 0.5 * (cubic(it->radius) + cubic((it - 1)->radius));
770  mean = pow(mean, 1. / 3.);
771  return mean;
772 }
773 
779 {
780  Mdouble rMin = getMaxRadius();
781  Mdouble rMax = getMinRadius();
782  return rMin / rMax;
783 }
784 
790 {
791  Mdouble SR = getSizeRatio();
792  if (SR > 100)
793  {
794  logger(WARN, "Size ratio > 100; Cutting the PSD to avoid inaccurate results");
795  cutoffCumulativeNumber(0.1, 0.9, 0.5);
796  }
797 }
798 
804 {
805  return particleSizeDistribution_[0].radius;
806 }
807 
813 {
814  return particleSizeDistribution_.back().radius;
815 }
816 
821 std::vector<PSD::RadiusAndProbability> PSD::getParticleSizeDistribution() const
822 {
824 }
825 
830 void PSD::setParticleSizeDistribution(std::vector<PSD::RadiusAndProbability> particleSizeDistribution)
831 {
832  particleSizeDistribution_ = particleSizeDistribution;
833 }
834 
841 {
842  int sum = 0;
843  for (auto& it: nParticlesPerClass_)
844  sum += it;
845  return sum;
846 }
847 
857 {
859  for (size_t im = 0; im < momenta_.size(); ++im)
860  {
861  // prevent summing up of moments for each time step of the simulation
862  momenta_[im] = 0;
863  for (auto& it : particleSizeDistribution_)
864  {
865  momenta_[im] += std::pow(it.radius, im) * it.probability;
866  }
867  }
868  // zeroth-moment equals particle number for a PROBABILITYDENSITY_NUMBER_DISTRIBUTION
871 }
872 
877 {
879  Mdouble mean = momenta_[1];
880  momenta_[5] += -5 * mean * momenta_[4] + 10 * mean * mean * momenta_[3]
881  - 10 * mean * mean * mean * momenta_[2] + 4 * mean * mean * mean * mean * mean;
882  momenta_[4] += -4 * mean * momenta_[3] + 6 * mean * mean * momenta_[2] - 3 * mean * mean * mean * mean;
883  momenta_[3] += -3 * mean * momenta_[2] + 2 * mean * mean * mean;
884  momenta_[2] += -mean * mean;
885 }
886 
891 {
893  Mdouble std = std::sqrt(momenta_[2]);
894  momenta_[3] /= std * std * std;
895  momenta_[4] /= std * std * std * std;
896  momenta_[5] /= std * std * std * std * std;
897 }
898 
903 std::array<Mdouble, 6> PSD::getMomenta() const
904 {
905  return momenta_;
906 }
907 
914 {
915  return l.probability < r.probability;
916 }
917 
924 bool operator<(const PSD::RadiusAndProbability& l, const Mdouble prob)
925 {
926  return l.probability < prob;
927 }
928 
934 {
935  return l.radius == r;
936 }
937 
942 std::istream& operator>>(std::istream& is, PSD::RadiusAndProbability& p)
943 {
944  is >> p.radius;
945  is >> p.probability;
946  return is;
947 }
948 
953 std::ostream& operator<<(std::ostream& os, PSD::RadiusAndProbability& p)
954 {
955  os << p.radius << ' ' << p.probability << ' ';
956  return os;
957 }
958 
PSD()
Constructor; sets everything to 0 or default.
Definition: PSD.cc:32
Class which stores radii and probabilities of a PSD. This class should be used as a vector
Definition: PSD.h:86
void computeStandardisedMomenta()
compute standardised momenta of the user defined PSD.
Definition: PSD.cc:890
MERCURY_DEPRECATED void setPSDFromVector(std::vector< RadiusAndProbability > psd, TYPE PSDType)
Deprecated version of reading in PSDs from a vector.
Definition: PSD.cc:352
std::vector< RadiusAndProbability > getParticleSizeDistribution() const
Get the PSD vector.
Definition: PSD.cc:821
~PSD()
Destructor; default destructor.
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here...
double Mdouble
Definition: GeneralDefine.h:34
void computeCentralMomenta()
compute central momenta of the user defined PSD.
Definition: PSD.cc:876
std::vector< int > nParticlesPerClass_
Definition: PSD.h:314
void setDistributionNormal(Mdouble mean, Mdouble standardDeviation, int numberOfBins)
create a PSD vector for a normal distribution.
Definition: PSD.cc:314
void convertCumulativeToCumulativeNumberDistribution(TYPE CDFType)
convert any other CDF to a CUMULATIVE_NUMBER_DISTRIBUTION.
Definition: PSD.cc:609
void convertProbabilityDensityNumberDistributionToProbabilityDensityVolumeDistribution()
convert a PROBABILITYDENSITY_NUMBER_DISTRIBUTION to a PROBABILITYDENSITY_VOLUME_DISTRIBUTION.
Definition: PSD.cc:584
const std::complex< Mdouble > i
Definition: ExtendedMath.h:51
Mdouble insertManuallyByVolume(Mdouble volume)
Draw sample radius manually per size class and check the volumeAllowed of each size class to insert t...
Definition: PSD.cc:137
TYPE
Enum class which stores the possible types of CDFs and PDFs. Particle size distributions can be repre...
Definition: PSD.h:71
void validateCumulativeDistribution()
Validates if a CDF starts with zero and adds up to unity.
Definition: PSD.cc:210
void read(const std::string &filename)
Reads .csv files into Mercury.
Definition: csvReader.cc:40
T cubic(const T val)
calculates the cube of a number
Definition: ExtendedMath.h:115
void setParticleSizeDistribution(std::vector< RadiusAndProbability >)
set the PSD by a suitable vector.
Definition: PSD.cc:830
#define MERCURY_DEPRECATED
Definition: GeneralDefine.h:37
void setHeader(bool headings)
Set the boolean hasHeader_.
Definition: csvReader.cc:109
PSD * copy() const
Creates a copy on the heap and returns a pointer.
Definition: PSD.cc:59
void cutoffAndSqueezeCumulative(Mdouble quantileMin, Mdouble quantileMax, Mdouble squeeze, Mdouble minPolydispersity=0.1)
cutoff the PSD at given quantiles and make it less polydisperse by squeezing it.
Definition: PSD.cc:702
std::array< Mdouble, 6 > momenta_
Definition: PSD.h:307
int getInsertedParticleNumber() const
Get the number of particles already inserted into the simulation.
Definition: PSD.cc:840
void validateProbabilityDensityDistribution()
Validates if the integral of the PDF equals to unity.
Definition: PSD.cc:251
std::vector< Mdouble > getFirstColumn(Mdouble scalingFactor)
Get first column of a .csv file and return it as a double.
Definition: csvReader.cc:128
Mdouble getRadiusByQuantile(Mdouble quantile) const
Calculate the quantile of the PSD.
Definition: PSD.cc:744
Mdouble getMinRadius() const
Get smallest radius of the PSD.
Definition: PSD.cc:803
std::istream & operator>>(std::istream &is, PSD::RadiusAndProbability &p)
Definition: PSD.cc:942
Mdouble drawSample()
Draw a sample radius from a CUMULATIVE_NUMBER_DISTRIBUTION.
Definition: PSD.cc:109
const Mdouble pi
Definition: ExtendedMath.h:45
void setDistributionUniform(Mdouble radMin, Mdouble radMax, int numberOfBins)
create a PSD vector for a uniform distribution.
Definition: PSD.cc:290
void cutoffCumulativeNumber(Mdouble quantileMin, Mdouble quantileMax, Mdouble minPolydispersity=0.1)
cutoff the PSD at given quantiles.
Definition: PSD.cc:668
std::vector< RadiusAndProbability > particleSizeDistribution_
Definition: PSD.h:301
std::vector< Mdouble > volumePerClass_
Definition: PSD.h:323
void printPSD() const
Prints radii and probabilities of the PSD vector.
Definition: PSD.cc:71
bool operator<(const PSD::RadiusAndProbability &l, const PSD::RadiusAndProbability &r)
Definition: PSD.cc:913
void convertProbabilityDensityToProbabilityDensityNumberDistribution(TYPE PDFType)
convert any PDF to a PROBABILITYDENSITY_NUMBER_DISTRIBUTION.
Definition: PSD.cc:522
void convertCumulativeToProbabilityDensity()
Converts a CDF to a PDF by derivation.
Definition: PSD.cc:502
Mdouble operator==(PSD::RadiusAndProbability l, const Mdouble r)
Definition: PSD.cc:933
Mdouble getVolumeDx(Mdouble x) const
Calculate a certain diameter (e.g. D10, D50, D90, etc.) from a percentile x of the volume based PSD...
Definition: PSD.cc:730
std::ostream & operator<<(std::ostream &os, PSD::RadiusAndProbability &p)
Definition: PSD.cc:953
Mdouble getVolumetricMeanRadius() const
get a volumetric mean radius of the PSD.
Definition: PSD.cc:762
std::vector< Mdouble > linspace(Mdouble a, Mdouble b, int N)
Definition: Helpers.cc:887
Mdouble getNumberDx(Mdouble x) const
Calculate a certain diameter (e.g. D10, D50, D90, etc.) from a percentile x of the number based PSD...
Definition: PSD.cc:720
Contains a vector with radii and probabilities of a user defined particle size distribution (PSD) ...
Definition: PSD.h:62
void setPSDFromCSV(const std::string &fileName, TYPE PSDType, bool headings=false, Mdouble unitScalingFactorRadii=1.0)
read in the PSD vector with probabilities and radii saved in a .csv file.
Definition: PSD.cc:418
Mdouble getSizeRatio() const
get the size ratio (width) of the PSD.
Definition: PSD.cc:778
void convertProbabilityDensityToCumulative()
Converts a PDF to a CDF by integration.
Definition: PSD.cc:487
T square(const T val)
squares a number
Definition: ExtendedMath.h:106
Mdouble getMaxRadius() const
Get largest radius of the PSD.
Definition: PSD.cc:812
std::vector< Mdouble > getSecondColumn(Mdouble scalingFactor)
Get second column of a .csv file and return it as a double.
Definition: csvReader.cc:145
Mdouble probability
Definition: PSD.h:90
std::array< Mdouble, 6 > getMomenta() const
get momenta of the user defined PSD.
Definition: PSD.cc:903
Enables reading of .csv files into MercuryDPM.
Definition: csvReader.h:45
void cutHighSizeRatio()
Check if the size ratio is too high and cut it.
Definition: PSD.cc:789
void computeRawMomenta()
compute raw momenta of the user defined PSD.
Definition: PSD.cc:856