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& p : momenta_)
35  p = 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 
72 {
73  if (particleSizeDistribution_.front().radius > 1e-1)
74  {
75  for (const auto p : particleSizeDistribution_)
76  {
77  std::cout << p.radius << "m\t" << p.probability * 100 << "\%\n";
78  }
79  }
80  else if (particleSizeDistribution_.front().radius > 1e-4)
81  {
82  for (const auto p : particleSizeDistribution_)
83  {
84  std::cout << p.radius * 1000 << "mm\t" << p.probability * 100 << "\%\n";
85  }
86  }
87  else if (particleSizeDistribution_.front().radius > 1e-7)
88  {
89  for (const auto p : particleSizeDistribution_)
90  {
91  std::cout << p.radius * 1e6 << "um\t" << p.probability * 100 << "\%\n";
92  }
93  }
94  else
95  {
96  for (const auto p : particleSizeDistribution_)
97  {
98  std::cout << p.radius * 1e9 << "nm\t" << p.probability * 100 << "\%\n";
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");
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(sum - 1 < 1e-6, "PDF is not valid: Integral of PDF is not equal to unity");
281 }
282 
295 void PSD::setPSDFromVector(std::vector<RadiusAndProbability> psdVector, TYPE PSDType)
296 {
297  particleSizeDistribution_ = psdVector;
298  switch (PSDType)
299  {
300  // default case is a cumulative number distribution (CND)
301  default:
303  break;
310  break;
317  break;
323  break;
328  break;
332  break;
337  break;
342  break;
343  }
344 }
345 
361 void PSD::setPSDFromCSV(const std::string& fileName, TYPE PSDType, bool headings, Mdouble unitScalingFactorRadii)
362 {
363  // read in csv file using the csvReader class
364  csvReader csv;
365  csv.setHeader(headings);
366  csv.read(fileName);
367  std::vector<Mdouble> radii = csv.getFirstColumn(unitScalingFactorRadii);
368  std::vector<Mdouble> probabilities = csv.getSecondColumn(1.0);
369  logger.assert_always(radii.size() == probabilities.size(), "The radii and probabilities vector have to be the "
370  "same size");
371  // combine radii and probabilities
372  for (int i = 0; i < radii.size(); ++i)
373  {
374  particleSizeDistribution_.push_back({radii[i], probabilities[i]});
375  }
376  logger.assert_always(!particleSizeDistribution_.empty(), "PSD cannot be empty");
377  switch (PSDType)
378  {
379  // default case is a CUMULATIVE_NUMBER_DISTRIBUTION
380  default:
382  break;
389  break;
395  break;
402  break;
406  break;
411  break;
416  break;
421  break;
422  }
423 }
424 
430 {
431  // add up probabilities
432  for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
433  {
434  it->probability = std::min(1.0, it->probability + (it - 1)->probability);
435  }
436  // check whether cumulative distribution is valid
438 }
439 
445 {
446  // subtract cumulative probabilities
447  Mdouble probabilityOld = 0, probabilityCDF;
448  for (auto& it : particleSizeDistribution_)
449  {
450  probabilityCDF = it.probability;
451  it.probability -= probabilityOld;
452  probabilityOld = probabilityCDF;
453  }
454  // check whether probability density distribution is valid
456 }
457 
465 {
466  Mdouble sum = 0;
467  switch (PDFType)
468  {
469  default:
470  logger(ERROR, "Wrong PDFType");
471  break;
473  for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
474  {
475  it->probability /= 0.5 * (it->radius + (it - 1)->radius);
476  // old conversion
477  //p.probability /= p.radius;
478  // sum up probabilities
479  sum += it->probability;
480  }
481  // normalize
482  for (auto& p : particleSizeDistribution_)
483  {
484  p.probability /= sum;
485  }
486  break;
488  for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
489  {
490  it->probability /=
491  1.0 / 3.0 * (square(it->radius) + it->radius * (it - 1)->radius + square((it - 1)->radius));
492  // old conversion
493  //p.probability /= square(p.radius);
494  // sum up probabilities
495  sum += it->probability;
496  }
497  // normalize
498  for (auto& p : particleSizeDistribution_)
499  {
500  p.probability /= sum;
501  }
502  break;
504  for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
505  {
506  it->probability /=
507  0.25 * (square(it->radius) + square((it - 1)->radius)) * (it->radius + (it - 1)->radius);
508  // old conversion
509  //p.probability /= cubic(p.radius);
510  // sum up probabilities
511  sum += it->probability;
512  }
513  // normalize
514  for (auto& p : particleSizeDistribution_)
515  {
516  p.probability /= sum;
517  }
518  break;
519  }
520 }
521 
527 {
528  Mdouble sum = 0;
529  for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
530  {
531  it->probability *=
532  0.25 * (square(it->radius) + square((it - 1)->radius)) * (it->radius + (it - 1)->radius);
533  // old conversion
534  //p.probability /= cubic(p.radius);
535  // sum up probabilities
536  sum += it->probability;
537  }
538  // normalize
539  for (auto& p : particleSizeDistribution_)
540  {
541  p.probability /= sum;
542  }
543 }
544 
552 {
553  Mdouble sum = 0;
554  switch (CDFType)
555  {
556  default:
557  logger(ERROR, "Wrong CDFType");
558  break;
560  for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
561  {
562  // add conversion here
563 
564  // sum up probabilities
565  sum += it->probability;
566  }
567  // normalize
568  for (auto& p : particleSizeDistribution_)
569  {
570  p.probability /= sum;
571  }
572  break;
574  for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
575  {
576  // add conversion here
577 
578  // sum up probabilities
579  sum += it->probability;
580  }
581  // normalize
582  for (auto& p : particleSizeDistribution_)
583  {
584  p.probability /= sum;
585  }
586  break;
588  for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
589  {
590  // add conversion here
591 
592  // sum up probabilities
593  sum += it->probability;
594  }
595  // normalize
596  for (auto& p : particleSizeDistribution_)
597  {
598  p.probability /= sum;
599  }
600  break;
601  }
602 }
603 
610 void PSD::cutoffCumulativeNumber(Mdouble percentileMin, Mdouble percentileMax, Mdouble minPolydispersity)
611 {
612  Mdouble radiusMin = getRadiusByPercentile(percentileMin);
613  Mdouble radiusMax = getRadiusByPercentile(percentileMax);
614  // to get a minimum polydispersity at the base
615  Mdouble radiusMinCut = std::min(radiusMin * (1 + minPolydispersity), radiusMax);
616  // cut off min
617  while (particleSizeDistribution_.front().radius <= radiusMinCut)
618  particleSizeDistribution_.erase(particleSizeDistribution_.begin());
619  particleSizeDistribution_.insert(particleSizeDistribution_.begin(),
620  particleSizeDistribution_[radiusMinCut, percentileMin]);
621  particleSizeDistribution_.insert(particleSizeDistribution_.begin(), particleSizeDistribution_[radiusMin, 0]);
622  // cut off max
623  while (particleSizeDistribution_.back().radius >= radiusMax) particleSizeDistribution_.pop_back();
624  particleSizeDistribution_.push_back(particleSizeDistribution_[radiusMax, percentileMax]);
625  particleSizeDistribution_.push_back(particleSizeDistribution_[radiusMax, 1]);
626 }
627 
637 void PSD::cutoffAndSqueezeCumulative(Mdouble percentileMin, Mdouble percentileMax, Mdouble squeeze, Mdouble
638 minPolydispersity)
639 {
640  Mdouble r50 = 0.5 * PSD::getDx(0.5);
641  // cut off
642  cutoffCumulativeNumber(percentileMin, percentileMax, minPolydispersity);
643  // squeeze psd
644  for (auto& p : particleSizeDistribution_)
645  {
646  p.radius = r50 + (p.radius - r50) * squeeze;
647  }
648 }
649 
656 {
657  return 2.0 * getRadiusByPercentile(x / 100);
658 }
659 
666 {
667  logger.assert_always(percentile <= 1 && percentile >= 0, "percentile is not between 0 and 1");
668  // find the percentile corresponding to the PSD
669  auto high = std::lower_bound(particleSizeDistribution_.begin(), particleSizeDistribution_.end(), percentile);
670  auto low = std::max(particleSizeDistribution_.begin(), high - 1);
671  if (high->probability == low->probability)
672  return high->radius;
673  else
674  return low->radius + (high->radius - low->radius) * (percentile - low->probability) /
675  (high->probability - low->probability);
676 }
677 
684 {
685  Mdouble mean = 0;
688  for (auto it = particleSizeDistribution_.begin() + 1; it != particleSizeDistribution_.end(); ++it)
689  mean += it->probability * 0.5 * (cubic(it->radius) + cubic((it - 1)->radius));
690  mean = pow(mean, 1. / 3.);
693  return mean;
694 }
695 
701 {
702  return particleSizeDistribution_[0].radius;
703 }
704 
710 {
711  return particleSizeDistribution_.back().radius;
712 }
713 
718 std::vector<PSD::RadiusAndProbability> PSD::getParticleSizeDistribution() const
719 {
721 }
722 
729 {
730  int sum = 0;
731  for (auto& it : nParticlesPerClass_)
732  sum += it;
733  return sum;
734 }
735 
745 {
747  for (size_t im = 0; im < momenta_.size(); ++im)
748  {
749  // prevent summing up of moments for each time step of the simulation
750  momenta_[im] = 0;
751  for (auto& it : particleSizeDistribution_)
752  {
753  momenta_[im] += std::pow(it.radius, im) * it.probability;
754  }
755  }
756  // zeroth-moment equals particle number for a PROBABILITYDENSITY_NUMBER_DISTRIBUTION
759 }
760 
765 {
767  Mdouble mean = momenta_[1];
768  momenta_[5] += -5 * mean * momenta_[4] + 10 * mean * mean * momenta_[3]
769  - 10 * mean * mean * mean * momenta_[2] + 4 * mean * mean * mean * mean * mean;
770  momenta_[4] += -4 * mean * momenta_[3] + 6 * mean * mean * momenta_[2] - 3 * mean * mean * mean * mean;
771  momenta_[3] += -3 * mean * momenta_[2] + 2 * mean * mean * mean;
772  momenta_[2] += -mean * mean;
773 }
774 
779 {
781  Mdouble std = std::sqrt(momenta_[2]);
782  momenta_[3] /= std * std * std;
783  momenta_[4] /= std * std * std * std;
784  momenta_[5] /= std * std * std * std * std;
785 }
786 
791 std::array<Mdouble, 6> PSD::getMomenta()
792 {
793  return momenta_;
794 }
795 
802 {
803  return l.probability < r.probability;
804 }
805 
812 bool operator<(const PSD::RadiusAndProbability& l, const Mdouble prob)
813 {
814  return l.probability < prob;
815 }
816 
822 {
823  return l.radius == r;
824 }
825 
830 std::istream& operator>>(std::istream& is, PSD::RadiusAndProbability& p)
831 {
832  is >> p.radius;
833  is >> p.probability;
834  return is;
835 }
836 
841 std::ostream& operator<<(std::ostream& os, PSD::RadiusAndProbability& p)
842 {
843  os << p.radius << ' ' << p.probability << ' ';
844  return os;
845 }
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:778
std::vector< RadiusAndProbability > getParticleSizeDistribution() const
Get the PSD vector.
Definition: PSD.cc:718
~PSD()
Destructor; default destructor.
MERCURY_DEPRECATED void setPSDFromVector(std::vector< RadiusAndProbability > psd, TYPE PSDType)
Deprecated version of reading in PSDs from a vector.
Definition: PSD.cc:295
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
double Mdouble
Definition: GeneralDefine.h:34
void computeCentralMomenta()
compute central momenta of the user defined PSD.
Definition: PSD.cc:764
std::vector< int > nParticlesPerClass_
Definition: PSD.h:284
void convertCumulativeToCumulativeNumberDistribution(TYPE CDFType)
convert any other CDF to a CUMULATIVE_NUMBER_DISTRIBUTION.
Definition: PSD.cc:551
void convertProbabilityDensityNumberDistributionToProbabilityDensityVolumeDistribution()
convert a PROBABILITYDENSITY_NUMBER_DISTRIBUTION to a PROBABILITYDENSITY_VOLUME_DISTRIBUTION.
Definition: PSD.cc:526
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
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:113
Mdouble getVolumetricMeanRadius()
get a volumetric mean radius of the PSD.
Definition: PSD.cc:683
Mdouble getMaxRadius()
Get largest radius of the PSD.
Definition: PSD.cc:709
Mdouble getDx(Mdouble x)
Calculate a certain diameter (e.g. D10, D50, D90, etc.) from a percentile of the PSD.
Definition: PSD.cc:655
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
std::array< Mdouble, 6 > momenta_
Definition: PSD.h:277
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
void cutoffAndSqueezeCumulative(Mdouble percentileMin, Mdouble percentileMax, Mdouble squeeze, Mdouble minPolydispersity=0.1)
cutoff the PSD at given percentiles and make it less polydisperse by squeezing it.
Definition: PSD.cc:637
std::istream & operator>>(std::istream &is, PSD::RadiusAndProbability &p)
Definition: PSD.cc:830
Mdouble drawSample()
Draw a sample radius from a CUMULATIVE_NUMBER_DISTRIBUTION.
Definition: PSD.cc:109
const Mdouble pi
Definition: ExtendedMath.h:45
Mdouble getMinRadius()
Get smallest radius of the PSD.
Definition: PSD.cc:700
std::vector< RadiusAndProbability > particleSizeDistribution_
Definition: PSD.h:272
void printPSD()
Prints radii and probabilities of the PSD vector.
Definition: PSD.cc:71
std::vector< Mdouble > volumePerClass_
Definition: PSD.h:293
bool operator<(const PSD::RadiusAndProbability &l, const PSD::RadiusAndProbability &r)
Definition: PSD.cc:801
void convertProbabilityDensityToProbabilityDensityNumberDistribution(TYPE PDFType)
convert any PDF to a PROBABILITYDENSITY_NUMBER_DISTRIBUTION.
Definition: PSD.cc:464
void cutoffCumulativeNumber(Mdouble percentileMin, Mdouble percentileMax, Mdouble minPolydispersity=0.1)
cutoff the PSD at given percentiles.
Definition: PSD.cc:610
void convertCumulativeToProbabilityDensity()
Converts a CDF to a PDF by derivation.
Definition: PSD.cc:444
Mdouble operator==(PSD::RadiusAndProbability l, const Mdouble r)
Definition: PSD.cc:821
int getInsertedParticleNumber()
Get the number of particles already inserted into the simulation.
Definition: PSD.cc:728
std::ostream & operator<<(std::ostream &os, PSD::RadiusAndProbability &p)
Definition: PSD.cc:841
std::array< Mdouble, 6 > getMomenta()
get momenta of the user defined PSD.
Definition: PSD.cc:791
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:361
Mdouble getRadiusByPercentile(Mdouble percentile)
Calculate the percentile of the PSD.
Definition: PSD.cc:665
void convertProbabilityDensityToCumulative()
Converts a PDF to a CDF by integration.
Definition: PSD.cc:429
T square(const T val)
squares a number
Definition: ExtendedMath.h:104
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
Enables reading of .csv files into MercuryDPM.
Definition: csvReader.h:45
void computeRawMomenta()
compute raw momenta of the user defined PSD.
Definition: PSD.cc:744