NAPISD
PAHdb website C++ backend
Loading...
Searching...
No Matches
main.cpp
1
2#include "../config.h"
3#include "Canvas.h"
4#include "FileInput.h"
5#include "FileOutput.h"
6#include "MinMax.h"
7#include "PAHEmissionModel.h"
8#include <exception>
9#include <new>
10#define MYSQLPP_SSQLS_NO_STATICS
11#include "PAHdb.h"
12#include "Parameters.h"
13#include "SpectralFitter.h"
14#include "ctime"
15#include <iostream>
16
18void new_handler() {
19
20 std::set_new_handler(nullptr);
21
22 time_t t = std::time(nullptr);
23
24 struct tm *tm_s = std::localtime(&t);
25
26 char str_t[32];
27
28 std::strftime(str_t, sizeof(str_t), "%Y%m%d-%X - EXCEPTION: ", tm_s);
29
30 std::cerr << str_t << "Memory allocation failed";
31
32 std::abort();
33}
34
42int main(const int argc, const char *argv[], char ** /* envp */) {
43
44 std::set_new_handler(new_handler);
45
46 Parameters parameters;
47
48 try {
49
50 parameters.parse(argc, argv);
51 } catch (const Exception &e) {
52
53 std::cerr << e.what() << std::endl;
54
55 return (1);
56 }
57
58 PAHdb pahdb;
59
60 try {
61
62 pahdb.connect(parameters.getDatabase(), parameters.getHost(),
63 parameters.getUsername(), parameters.getPassword(),
64 parameters.getPort(), parameters.isCompress(),
65 parameters.getTimeout(), parameters.getSocket());
66
67 pahdb.setTable(parameters.getTable());
68
69 pahdb.setProgress(10.0);
70 } catch (const Exception &e) {
71
72 std::cerr << e.what() << std::endl;
73
74 return (1);
75 }
76
77 std::vector<std::vector<std::pair<double, double>>> transitions;
78
79 std::vector<std::vector<double>> spectra;
80
81 std::vector<double> grid;
82
83 std::vector<double> temperatures;
84
85 std::tuple<std::vector<double>, std::vector<double>, std::vector<double>>
86 inputdata;
87
88 PAHEmissionModel pahemissionmodel;
89
90 SpectralFitter spectralfitter;
91
92 try {
93
94 if (parameters.getTool() == Parameters::Arg::CompareExperimentWithTheory) {
95
96 transitions = pahdb.getExperimentalAndTheoreticalTransitionsFromId(
97 parameters.getIds().at(0));
98 } else {
99
100 if (parameters.getTool() == Parameters::Arg::SpectralFit) {
101
102 FileInput fileinput(parameters.getInputFilename());
103
104 inputdata = fileinput.readFileFromDisk();
105
106 pahemissionmodel.setGrid(std::get<0>(inputdata));
107
108 std::array<double, 2> inputlimits;
109
110 inputlimits[0] = *max_element(std::get<0>(inputdata).begin(),
111 std::get<0>(inputdata).end());
112
113 inputlimits[1] = *min_element(std::get<0>(inputdata).begin(),
114 std::get<0>(inputdata).end());
115
116 if (parameters.getUnits() == Parameters::Arg::Wavelength) {
117
118 PAHEmissionModel::convertFromFrequencyToWavelength(inputlimits);
119 }
120
121 parameters.setPlotLimits(inputlimits);
122
123 spectralfitter.setSpectrum(std::get<1>(inputdata));
124
125 spectralfitter.setUncertainties(std::get<2>(inputdata));
126
127 } else {
128
129 pahemissionmodel.makeGrid(parameters.getPlotLimitsAsFrequencies().at(0),
130 parameters.getPlotLimitsAsFrequencies().at(1),
131 2.0);
132 }
133
134 grid = pahemissionmodel.getGrid();
135
136 if (parameters.getTool() != Parameters::Arg::TemperatureStack) {
137
138 pahemissionmodel.setTransitions(
139 pahdb.getTransitionsFromIds(parameters.getIds()));
140 } else {
141
142 transitions = pahdb.getTransitionsFromId(parameters.getIds().at(0));
143
144 transitions = std::vector<std::vector<std::pair<double, double>>>(
145 parameters.getTemperaturesInKelvin().size(), transitions.at(0));
146
147 pahemissionmodel.setTransitions(transitions);
148 }
149
150 if (parameters.getTable() == PAHdb::Database::Experiment &&
151 (parameters.getModel() == Parameters::Arg::Temperature ||
152 parameters.getModel() == Parameters::Arg::Cascade)) {
153 std::vector<sql_properties> prop = pahdb.getPropertiesByIDs(
154 parameters.getIds()); // perhaps other name then 'prop'?
155
156 std::vector<int> charges;
157 charges.reserve(parameters.getIds().size());
158 std::vector<int> carbons;
159 carbons.reserve(parameters.getIds().size());
160 for (const auto &properties : prop) {
161
162 charges.push_back(properties.charge);
163 carbons.push_back(properties.n_c);
164 }
165
166 pahemissionmodel.useApproximate(charges, carbons);
167 }
168
169 switch (parameters.getModel()) {
170
171 case Parameters::Arg::FixedTemperature:
172
173 if (parameters.getTool() == Parameters::Arg::TemperatureStack) {
174
175 pahemissionmodel.applyBlackbodyWithTemperatureForEach(
176 parameters.getTemperaturesInKelvin());
177 } else {
178
179 pahemissionmodel.applyBlackbodyWithTemperature(
180 parameters.getTemperatureInKelvin());
181 }
182
183 break;
184
185 case Parameters::Arg::Temperature:
186
187 pahemissionmodel.applyTemperatureWithEnergy(parameters.getEnergyInCGS(),
188 temperatures);
189
190 break;
191
192 case Parameters::Arg::Cascade:
193
194 pahemissionmodel.applyCascadeWithEnergy(parameters.getEnergyInCGS(),
195 temperatures);
196
197 break;
198
199 case Parameters::Arg::ZeroKelvin:
200 default:
201
202 break;
203 };
204
205 if (parameters.getShift() != 0.0) {
206
207 pahemissionmodel.shiftTransitions(parameters.getShift());
208 }
209
210 switch (parameters.getProfile()) {
211
212 case Parameters::Arg::Gaussian:
213
214 pahemissionmodel.getSpectraAndConvolveWithGaussianOfFHWM(
215 spectra, parameters.getFWHM());
216
217 break;
218 case Parameters::Arg::Drude:
219
220 pahemissionmodel.getSpectraAndConvolveWithDrudeOfFHWM(
221 spectra, parameters.getFWHM());
222
223 break;
224 case Parameters::Arg::Lorentzian:
225 default:
226
227 pahemissionmodel.getSpectraAndConvolveWithLorentianOfFHWM(
228 spectra, parameters.getFWHM());
229 break;
230 };
231 }
232
233 } catch (const Exception &e) {
234
235 std::cerr << e.what() << std::endl;
236
237 try {
238
239 pahdb.setError(e.what());
240 } catch (const Exception &e) {
241
242 std::cerr << e.what() << std::endl;
243 }
244
245 return (1);
246 }
247
248 try {
249
250 pahdb.setProgress(60.0);
251 } catch (const Exception &e) {
252
253 std::cerr << e.what() << std::endl;
254
255 return (1);
256 }
257
258 std::vector<double> small;
259
260 std::vector<double> large;
261
262 std::vector<double> anion;
263
264 std::vector<double> neutral;
265
266 std::vector<double> cation;
267
268 if (parameters.getTool() == Parameters::Arg::SpectralFit) {
269
270 spectralfitter.setPool(spectra);
271
272 spectralfitter.fitSpectrum();
273
274 std::vector<unsigned int> indices = spectralfitter.getNonZeroIndices();
275
276 std::vector<int> ids;
277
278 for (const auto index : indices) {
279
280 ids.push_back(parameters.getIds().at(index));
281 }
282
283 std::vector<sql_properties> prop =
284 pahdb.getPropertiesByIDs(ids); // perhaps other name then 'prop'?
285
286 small = std::vector<double>(grid.size(), 0.0);
287
288 large = std::vector<double>(grid.size(), 0.0);
289
290 anion = std::vector<double>(grid.size(), 0.0);
291
292 neutral = std::vector<double>(grid.size(), 0.0);
293
294 cation = std::vector<double>(grid.size(), 0.0);
295
296 unsigned int j = 0;
297
298 for (const auto &properties : prop) {
299
300 unsigned int k = 0;
301
302 if (properties.n_c > 50) {
303
304 for (auto &l : large) {
305
306 l += spectralfitter.getWeights().at(indices.at(j)) *
307 spectra.at(indices.at(j)).at(k++);
308 }
309 } else {
310
311 for (auto &s : small) {
312
313 s += spectralfitter.getWeights().at(indices.at(j)) *
314 spectra.at(indices.at(j)).at(k++);
315 }
316 }
317
318 k = 0;
319
320 if (properties.charge < 0) {
321
322 for (auto &a : anion) {
323
324 a += spectralfitter.getWeights().at(indices.at(j)) *
325 spectra.at(indices.at(j)).at(k++);
326 }
327 } else if (properties.charge == 0) {
328
329 for (auto &n : neutral) {
330
331 n += spectralfitter.getWeights().at(indices.at(j)) *
332 spectra.at(indices.at(j)).at(k++);
333 }
334 } else {
335
336 for (auto &c : cation) {
337
338 c += spectralfitter.getWeights().at(indices.at(j)) *
339 spectra.at(indices.at(j)).at(k++);
340 }
341 }
342
343 ++j;
344 }
345 }
346
347 try {
348
349 pahdb.setProgress(80.0);
350 } catch (const Exception &e) {
351
352 std::cerr << e.what() << std::endl;
353
354 return (1);
355 }
356
357 Canvas canvas;
358
359 canvas.setSize(parameters.getPlotSize());
360
361 canvas.setColor("ffffff");
362
363 canvas.setDefaultCharacterHeight(parameters.getPlotSize().at(1) * 2.0 /
364 325.0);
365
366 Panels panels;
367
368 Plot plot;
369
370 plot.setXLimits(parameters.getPlotLimits());
371
372 plot.setDrawHorizontalFineGrid();
373
374 plot.setDrawVerticalFineGrid();
375
376 plot.getHorizontalGrid().setColor("dedede");
377
378 plot.getHorizontalGrid().setLineStyle(Plot::Style::Dashed);
379
380 plot.getVerticalGrid() = plot.getHorizontalGrid();
381
382 plot.getYAxis().at(0).setWriteLabelsParallel();
383
384 plot.getYAxis().at(0).setStyle(Axis::Style::WritePowerInFrame);
385
386 std::string xtitle;
387
388 std::string xtitle2;
389
390 std::string ytitle;
391
392 switch (parameters.getUnits()) {
393
394 case Parameters::Arg::Wavelength:
395
396 xtitle = "#fiwavelength [micron]";
397
398 xtitle2 = "#fifrequency [cm#u-1#d]";
399
400 if (parameters.getTool() != Parameters::Arg::CompareExperimentWithTheory) {
401
402 PAHEmissionModel::convertFromFrequencyToWavelength(grid);
403 } else {
404
405 PAHEmissionModel::convertFromFrequencyToWavelength(transitions);
406 }
407
408 break;
409 case Parameters::Arg::Frequency:
410 default:
411
412 xtitle = "#fifrequency [cm#u-1#d]";
413
414 xtitle2 = "#fiwavelength [micron]";
415 }
416
417 if (parameters.getTool() == Parameters::Arg::CompareExperimentWithTheory) {
418
419 ytitle = "#fiintegrated cross-section [#fnx#fi10#u5#d cm/mol]";
420 } else if (parameters.getModel() == Parameters::Arg::ZeroKelvin) {
421
422 ytitle = "#ficross-section [#fnx#fi10#u5#d cm#u2#d/mol]";
423 } else if (parameters.getModel() == Parameters::Arg::FixedTemperature ||
424 parameters.getModel() == Parameters::Arg::Temperature) {
425
426 ytitle = "#fispectral radiance [erg/s/cm#u-1#d/PAH]";
427 } else if (parameters.getModel() == Parameters::Arg::Cascade) {
428
429 ytitle = "#firadiant energy [erg/cm#u-1#d/PAH]";
430 }
431
432 std::vector<Plot> plots;
433
434 plots.reserve(2);
435
436 Curve curve;
437
438 curve.setFill();
439
440 Line line;
441
442 Point point;
443
444 point.setSymbol(256);
445
446 Text text;
447
448 text.setJustification(Text::RightJustification);
449
450 text.setColor("0000ff");
451
452 std::vector<std::string> formulae;
453
454 formulae.reserve(2);
455
456 std::vector<PAHGeometry> geometries;
457
458 geometries.reserve(2);
459
460 std::vector<double> sum;
461
462 std::ostringstream ostr;
463
464 std::array<double, 2> margins;
465
466 std::array<double, 2> newmargins;
467
468 int i = 0, j = 0;
469
470 std::array<int, 2> canvassize;
471
472 std::vector<bool> available;
473
474 available.reserve(3);
475
476 int npanels;
477
478 std::array<const char *, 3> annotations = {"Laboratory", "Harmonic Theory",
479 "Anharmonic Theory"};
480
481 switch (parameters.getTool()) {
482
483 case Parameters::Arg::Stack:
484 case Parameters::Arg::TemperatureStack:
485
486 canvassize[0] = canvas.getSize().at(0);
487
488 canvassize[1] = canvas.getSize().at(1) +
489 (spectra.size() - 1) * parameters.getPlotSize().at(1) *
490 (plot.getYMargins().at(1) - plot.getYMargins().at(0));
491
492 canvas.setSize(canvassize);
493
494 margins[0] = plot.getYMargins().at(0) * parameters.getPlotSize().at(1) /
495 canvas.getSize().at(1);
496
497 margins[1] = 1.0 - ((1.0 - plot.getYMargins().at(1)) *
498 parameters.getPlotSize().at(1)) /
499 canvas.getSize().at(1);
500
501 panels.setYMargins(margins);
502
503 panels.setLayout(1, spectra.size());
504
505 geometries = pahdb.getGeometriesFromIds(parameters.getIds());
506
507 std::for_each(geometries.begin(), geometries.end(),
508 std::mem_fn(&PAHGeometry::diagonalize));
509
510 for (const auto &spectrum : spectra) {
511
512 plot.setAdvance(false);
513
514 // Draw panel grid
515
516 plot.setDrawHorizontalFineGrid();
517
518 plot.setDrawVerticalFineGrid();
519
520 plot.getXAxis().at(0).setDrawConventionalAxis(false);
521
522 plot.getXAxis().at(0).setDrawUnconventionalAxis(false);
523
524 plot.getXAxis().at(0).setDrawConventionalLabels(false);
525
526 plot.getYAxis().at(0).setDrawConventionalAxis(false);
527
528 plot.getYAxis().at(0).setDrawUnconventionalAxis(false);
529
530 plot.getYAxis().at(0).setDrawConventionalLabels(false);
531
532 plot.setYLimits(
533 MinMax::min_max(spectrum, MinMax::Nice | MinMax::MaxExtraRoom));
534
535 plots.push_back(plot);
536
537 plot.clear();
538
539 plot.setDrawHorizontalFineGrid(false);
540
541 plot.setDrawVerticalFineGrid(false);
542
543 // Draw panel structure
544
545 if (parameters.getTool() == Parameters::Arg::Stack) {
546
547 line.setColor("cccccc");
548
549 line.setLineWidth(5);
550
551 j = 0;
552
553 for (const auto &bonds : geometries.at(i).getBonds()) {
554
555 for (const auto &bond : bonds) {
556
557 line.setStartCoordinates(
558 (plot.getXLimits().at(1) + plot.getXLimits().at(0)) / 2.0 +
559 geometries.at(i)[j].getX() *
560 (plot.getXLimits().at(1) - plot.getXLimits().at(0)) /
561 (28.0 * parameters.getPlotSize().at(0) /
562 parameters.getPlotSize().at(1)),
563 (plot.getYLimits().at(1) + plot.getYLimits().at(0)) / 2.0 +
564 geometries.at(i)[j].getY() *
565 (plot.getYLimits().at(1) - plot.getYLimits().at(0)) /
566 28.0,
567 geometries.at(i)[j].getZ());
568
569 line.setEndCoordinates(
570 (plot.getXLimits().at(1) + plot.getXLimits().at(0)) / 2.0 +
571 geometries.at(i)[bond].getX() *
572 (plot.getXLimits().at(1) - plot.getXLimits().at(0)) /
573 (28.0 * parameters.getPlotSize().at(0) /
574 parameters.getPlotSize().at(1)),
575 (plot.getYLimits().at(1) + plot.getYLimits().at(0)) / 2.0 +
576 geometries.at(i)[bond].getY() *
577 (plot.getYLimits().at(1) - plot.getYLimits().at(0)) /
578 28.0,
579 geometries.at(i)[j].getZ());
580
581 plot.add(line);
582 }
583
584 ++j;
585 }
586
587 for (const auto &atom : geometries.at(i)) {
588
589 point.setCoordinates(
590 (plot.getXLimits().at(1) + plot.getXLimits().at(0)) / 2.0 +
591 atom.getX() *
592 (plot.getXLimits().at(1) - plot.getXLimits().at(0)) /
593 (28.0 * parameters.getPlotSize().at(0) /
594 parameters.getPlotSize().at(1)),
595 (plot.getYLimits().at(1) + plot.getYLimits().at(0)) / 2.0 +
596 atom.getY() *
597 (plot.getYLimits().at(1) - plot.getYLimits().at(0)) /
598 28.0,
599 atom.getZ());
600
601 point.setColor(atom.getCPKColor());
602
603 point.setSize(atom.getSize() / 70.0);
604
605 plot.add(point);
606 }
607
608 plots.push_back(plot);
609
610 plot.clear();
611 }
612
613 // Write panel formula
614
615 formulae = pahdb.getFormulaeFromIds(parameters.getIds());
616
617 if (parameters.getTool() != Parameters::Arg::TemperatureStack) {
618
619 text.setText(Text::formatChemicalFormula(formulae.at(i)));
620
621 text.setCoordinates(0.95, 0.85, 0.0, Text::CoordinateSystem::NORMAL);
622
623 plot.add(text);
624
625 plots.push_back(plot);
626
627 plot.clear();
628 }
629
630 // Write panel temperature
631
632 if ((parameters.getModel() != Parameters::Arg::ZeroKelvin &&
633 parameters.getModel() != Parameters::Arg::FixedTemperature) ||
634 parameters.getTool() == Parameters::Arg::TemperatureStack) {
635
636 ostr.str("");
637
638 ostr << static_cast<int>(
639 parameters.getTool() == Parameters::Arg::TemperatureStack
640 ? parameters.getTemperaturesInKelvin().at(i)
641 : temperatures.at(i))
642 << " K";
643
644 text.setText(ostr.str());
645
646 text.setCoordinates(0.95, 0.75, 0.0, Text::CoordinateSystem::NORMAL);
647
648 plot.add(text);
649
650 plots.push_back(plot);
651
652 plot.clear();
653 }
654
655 plot.getXAxis().at(0).setDrawConventionalAxis();
656
657 plot.getXAxis().at(0).setDrawUnconventionalAxis();
658
659 plot.getYAxis().at(0).setDrawConventionalAxis();
660
661 plot.getYAxis().at(0).setDrawUnconventionalAxis();
662
663 plot.getYAxis().at(0).setDrawConventionalLabels();
664
665 // Draw panel spectrum
666
667 curve.setXAndY(grid, spectra.at(i++));
668
669 plot.add(curve);
670
671 plot.setAdvance(true);
672
673 plots.push_back(plot);
674
675 plot.clear();
676 }
677
678 plots.front().getXAxis().at(0).setDrawConventionalLabels();
679
680 plots.front().getXAxis().at(0).setTitle(xtitle);
681
682 plots.back().getXAxis().at(0).setDrawUnconventionalAxis(false);
683
684 plots.back().getXAxis().emplace_back();
685
686 plots.back().getXAxis().at(1).setDrawConventionalLabels(false);
687
688 // plots.back().getXAxis().at(1).setReciprocalTickFinder();
689
690 plots.back().getXAxis().at(1).setReciprocalLabelFormatter();
691
692 plots.back().getXAxis().at(1).setDrawConventionalAxis(false);
693
694 plots.back().getXAxis().at(1).setDrawUnconventionalLabels();
695
696 plots.back().getXAxis().at(1).setTitle(xtitle2);
697
698 panels.add(plots);
699
700 canvas.add(panels);
701
702 text.setColor("000000");
703
704 text.setText(ytitle);
705
706 text.setJustification(Text::CenterJustification);
707
708 text.setAngle(90.0);
709
710 text.setSize(1.0);
711
712 text.setCoordinates(
713 0.0723,
714 panels.getYMargins().at(0) +
715 (panels.getYMargins().at(1) - panels.getYMargins().at(0)) / 2.0,
716 0.0);
717
718 canvas.add(text);
719
720 break;
721 case Parameters::Arg::CompareExperimentWithTheory:
722
723 for (const auto &t : transitions) {
724
725 available.emplace_back(t.size());
726 }
727
728 npanels = std::count(available.begin(), available.end(), true);
729
730 canvassize[0] = canvas.getSize().at(0);
731
732 canvassize[1] = canvas.getSize().at(1) +
733 (npanels - 1) * parameters.getPlotSize().at(1) *
734 (plot.getYMargins().at(1) - plot.getYMargins().at(0));
735
736 canvas.setSize(canvassize);
737
738 margins[0] = plot.getYMargins().at(0) * parameters.getPlotSize().at(1) /
739 canvas.getSize().at(1);
740
741 margins[1] = 1.0 - ((1.0 - plot.getYMargins().at(1)) *
742 parameters.getPlotSize().at(1)) /
743 canvas.getSize().at(1);
744
745 // Draw panel grid
746
747 plot.setAdvance(false);
748
749 plot.getXAxis().at(0).setDrawConventionalAxis(false);
750
751 plot.getXAxis().at(0).setDrawUnconventionalAxis(false);
752
753 plot.getXAxis().at(0).setDrawConventionalLabels(false);
754
755 plot.getYAxis().at(0).setDrawConventionalAxis(false);
756
757 plot.getYAxis().at(0).setDrawUnconventionalAxis(false);
758
759 plot.getYAxis().at(0).setDrawConventionalLabels(false);
760
761 plot.setYLimits(
762 MinMax::min_max(transitions, MinMax::Nice | MinMax::MaxExtraRoom));
763
764 newmargins[0] = margins[0];
765
766 newmargins[1] = margins[0] + (margins[1] - margins[0]) / npanels;
767
768 plot.setYMargins(newmargins);
769
770 plots.push_back(plot);
771
772 plot.clear();
773
774 newmargins[0] = newmargins[1];
775
776 newmargins[1] = newmargins[1] + (margins[1] - margins[0]) / npanels;
777
778 plot.setYMargins(newmargins);
779
780 plots.push_back(plot);
781
782 plot.clear();
783
784 if (npanels > 2) {
785
786 newmargins[0] = newmargins[1];
787
788 newmargins[1] = newmargins[1] + (margins[1] - margins[0]) / npanels;
789
790 plot.setYMargins(newmargins);
791
792 plots.push_back(plot);
793
794 plot.clear();
795 }
796
797 // Draw canvas geometry
798
799 geometries = pahdb.getGeometriesFromIds(parameters.getIds());
800
801 std::for_each(geometries.begin(), geometries.end(),
802 std::mem_fn(&PAHGeometry::diagonalize));
803
804 plot.setDrawHorizontalFineGrid(false);
805
806 plot.setDrawVerticalFineGrid(false);
807
808 plot.setYMargins(margins);
809
810 line.setColor("cccccc");
811
812 line.setLineWidth(5);
813
814 for (const auto &bonds : geometries.at(0).getBonds()) {
815
816 for (const auto bond : bonds) {
817
818 line.setStartCoordinates(
819 (plot.getXLimits().at(1) + plot.getXLimits().at(0)) / 2.0 +
820 geometries.at(0)[i].getX() *
821 (plot.getXLimits().at(1) - plot.getXLimits().at(0)) / 28.0,
822 (plot.getYLimits().at(1) + plot.getYLimits().at(0)) / 2.0 +
823 canvas.getAspectRatio() * geometries.at(0)[i].getY() *
824 (plot.getYLimits().at(1) - plot.getYLimits().at(0)) / 28.0,
825 geometries.at(0)[i].getZ());
826
827 line.setEndCoordinates(
828 (plot.getXLimits().at(1) + plot.getXLimits().at(0)) / 2.0 +
829 geometries.at(0)[bond].getX() *
830 (plot.getXLimits().at(1) - plot.getXLimits().at(0)) / 28.0,
831 (plot.getYLimits().at(1) + plot.getYLimits().at(0)) / 2.0 +
832 canvas.getAspectRatio() * geometries.at(0)[bond].getY() *
833 (plot.getYLimits().at(1) - plot.getYLimits().at(0)) / 28.0,
834 geometries.at(0)[bond].getZ());
835
836 plot.add(line);
837 }
838
839 ++i;
840 }
841
842 i = 0;
843
844 for (const auto &atom : geometries.at(0)) {
845
846 point.setCoordinates(
847 (plot.getXLimits().at(1) + plot.getXLimits().at(0)) / 2.0 +
848 atom.getX() *
849 (plot.getXLimits().at(1) - plot.getXLimits().at(0)) / 28.0,
850 (plot.getYLimits().at(1) + plot.getYLimits().at(0)) / 2.0 +
851 canvas.getAspectRatio() * atom.getY() *
852 (plot.getYLimits().at(1) - plot.getYLimits().at(0)) / 28.0,
853 atom.getZ());
854
855 point.setColor(atom.getCPKColor());
856
857 point.setSize(atom.getSize() / 40.0);
858
859 plot.add(point);
860 }
861
862 // Draw canvas formula
863
864 formulae = pahdb.getFormulaeFromIds(parameters.getIds());
865
866 text.setColor("0000ff");
867
868 text.setSize(2);
869
870 text.setJustification(Text::RightJustification);
871
872 text.setText(Text::formatChemicalFormula(formulae.at(0)));
873
874 text.setCoordinates(0.95, 0.95, 0.0, Text::CoordinateSystem::NORMAL);
875
876 plot.add(text);
877
878 plots.push_back(plot);
879
880 plot.clear();
881
882 plot.getXAxis().at(0).setTitle(xtitle);
883
884 plot.getXAxis().at(0).setDrawConventionalAxis();
885
886 if (npanels > 2) {
887
888 plot.getXAxis().emplace_back();
889
890 plot.getXAxis().at(1).setDrawConventionalLabels(false);
891 }
892
893 plot.getXAxis().at(0).setDrawConventionalLabels();
894
895 plot.getYAxis().at(0).setDrawConventionalAxis();
896
897 plot.getYAxis().at(0).setDrawUnconventionalAxis();
898
899 plot.getYAxis().at(0).setDrawConventionalLabels();
900
901 newmargins[0] = margins[0];
902
903 newmargins[1] = margins[0] + (margins[1] - margins[0]) / npanels;
904
905 plot.setYMargins(newmargins);
906
907 line.setLineWidth(4);
908
909 line.setColor("ff0000");
910
911 if (!available.at(i)) {
912
913 ++i;
914 }
915
916 // Draw panel annotation
917
918 text.setSize(3);
919
920 text.setJustification(Text::CenterJustification);
921
922 text.setText(annotations.at(i));
923
924 text.setCoordinates(0.5, 0.75, 0.0, Text::CoordinateSystem::NORMAL);
925
926 plot.add(text);
927
928 // Draw panel transitions
929
930 for (const auto &transition : transitions.at(i++)) {
931
932 line.setStartCoordinates(transition.first, 0.0, 0.0);
933
934 line.setEndCoordinates(transition.first, transition.second, 0.0);
935
936 plot.add(line);
937 }
938
939 plots.push_back(plot);
940
941 plot.clear();
942
943 newmargins[0] = newmargins[1];
944
945 newmargins[1] = newmargins[1] + (margins[1] - margins[0]) / npanels;
946
947 plot.setYMargins(newmargins);
948
949 if (npanels > 2) {
950
951 plot.getXAxis().at(0).setTitle("");
952
953 plot.getXAxis().at(0).setDrawConventionalLabels(false);
954 }
955
956 if (npanels < 3) {
957
958 plot.getXAxis().emplace_back();
959
960 // plot.getXAxis().at(1).setReciprocalTickFinder();
961
962 plot.getXAxis().at(1).setReciprocalLabelFormatter();
963
964 plot.getXAxis().at(0).setDrawUnconventionalAxis(false);
965
966 plot.getXAxis().at(1).setDrawConventionalAxis(false);
967
968 plot.getXAxis().at(1).setDrawConventionalLabels(false);
969
970 plot.getXAxis().at(1).setDrawUnconventionalAxis();
971
972 plot.getXAxis().at(1).setDrawUnconventionalLabels();
973
974 plot.getXAxis().at(1).setTitle(xtitle2);
975
976 plot.getXAxis().at(0).setTitle("");
977
978 plot.getXAxis().at(0).setDrawZeroLine();
979
980 plot.getXAxis().at(0).setDrawConventionalLabels(false);
981 }
982
983 if (!available.at(i)) {
984
985 ++i;
986 }
987
988 // Draw panel annotation
989
990 text.setText(annotations.at(i));
991
992 text.setCoordinates(0.5, 0.75, 0.0, Text::CoordinateSystem::NORMAL);
993
994 plot.add(text);
995
996 // Draw panel transitions
997
998 for (const auto &transition : transitions.at(i++)) {
999
1000 line.setStartCoordinates(transition.first, 0.0, 0.0);
1001
1002 line.setEndCoordinates(transition.first, transition.second, 0.0);
1003
1004 plot.add(line);
1005 }
1006
1007 plots.push_back(plot);
1008
1009 plot.clear();
1010
1011 if (npanels > 2) {
1012
1013 newmargins[0] = newmargins[1];
1014
1015 newmargins[1] = newmargins[1] + (margins[1] - margins[0]) / npanels;
1016
1017 plot.setYMargins(newmargins);
1018
1019 // plot.getXAxis().at(1).setReciprocalTickFinder();
1020
1021 plot.getXAxis().at(1).setReciprocalLabelFormatter();
1022
1023 plot.getXAxis().at(0).setDrawUnconventionalAxis(false);
1024
1025 plot.getXAxis().at(1).setDrawConventionalAxis(false);
1026
1027 plot.getXAxis().at(1).setDrawConventionalLabels(false);
1028
1029 plot.getXAxis().at(1).setDrawUnconventionalAxis();
1030
1031 plot.getXAxis().at(1).setDrawUnconventionalLabels();
1032
1033 plot.getXAxis().at(1).setTitle(xtitle2);
1034
1035 plot.getXAxis().at(0).setTitle("");
1036
1037 plot.getXAxis().at(0).setDrawZeroLine();
1038
1039 plot.getXAxis().at(0).setDrawConventionalLabels(false);
1040
1041 // Draw panel annotation
1042
1043 text.setText(annotations.at(i));
1044
1045 text.setCoordinates(0.5, 0.75, 0.0, Text::CoordinateSystem::NORMAL);
1046
1047 plot.add(text);
1048
1049 // Draw panel transitions
1050
1051 for (const auto &transition : transitions.at(i)) {
1052
1053 line.setStartCoordinates(transition.first, 0.0, 0.0);
1054
1055 line.setEndCoordinates(transition.first, transition.second, 0.0);
1056
1057 plot.add(line);
1058 }
1059
1060 plots.push_back(plot);
1061 }
1062
1063 canvas.add(plots);
1064
1065 text.setColor("000000");
1066
1067 text.setText(ytitle);
1068
1069 text.setJustification(Text::CenterJustification);
1070
1071 text.setAngle(90.0);
1072
1073 text.setSize(1.0);
1074
1075 text.setCoordinates(0.0723,
1076 margins[0] + (margins[1] - margins[0]) / npanels, 0.0);
1077
1078 canvas.add(text);
1079
1080 break;
1081 case Parameters::Arg::SpectralFit:
1082
1083 curve.setSymbol(9);
1084
1085 curve.setXAndY(grid, std::get<1>(inputdata));
1086
1087 if (std::get<2>(inputdata).size() == std::get<1>(inputdata).size()) {
1088
1089 curve.setYErr(std::get<2>(inputdata));
1090 }
1091
1092 plot.add(curve);
1093
1094 curve.clear();
1095
1096 curve.setSymbol(0);
1097
1098 curve.setFill(false);
1099
1100 curve.setColor("dedede");
1101
1102 sum = std::vector<double>(spectra.at(0).size());
1103
1104 for (auto &spectrum : spectra) {
1105
1106 int j = 0;
1107
1108 if (spectralfitter.getWeights().at(i) == 0.0) {
1109
1110 ++i;
1111
1112 continue;
1113 }
1114
1115 for (auto &s : spectrum) {
1116
1117 s *= spectralfitter.getWeights().at(i);
1118
1119 sum.at(j++) += s;
1120 }
1121
1122 curve.setXAndY(grid, spectrum);
1123
1124 plot.add(curve);
1125
1126 ++i;
1127 }
1128
1129 curve.setXAndY(grid, sum);
1130
1131 curve.setColor("ff0000");
1132
1133 curve.setLineWidth(2);
1134
1135 plot.add(curve);
1136
1137 plot.setYLimits(MinMax::min_max(sum, MinMax::Nice));
1138
1139 text.setJustification(Text::LeftJustification);
1140
1141 text.setColor("ff0000");
1142
1143 text.setText("total fit");
1144
1145 text.setCoordinates(0.85, 0.90, 0.0, Text::CoordinateSystem::NORMAL);
1146
1147 plot.add(text);
1148
1149 text.setColor("00ff00");
1150
1151 text.setText("N#dC#u>50");
1152
1153 text.setCoordinates(0.85, 0.83, 0.0, Text::CoordinateSystem::NORMAL);
1154
1155 plot.add(text);
1156
1157 text.setColor("0000ff");
1158
1159 text.setText("N#dC#u<50");
1160
1161 text.setCoordinates(0.85, 0.76, 0.0, Text::CoordinateSystem::NORMAL);
1162
1163 plot.add(text);
1164
1165 text.setColor("fff00f");
1166
1167 text.setText("anion");
1168
1169 text.setCoordinates(0.85, 0.69, 0.0, Text::CoordinateSystem::NORMAL);
1170
1171 plot.add(text);
1172
1173 text.setColor("ff00ff");
1174
1175 text.setText("neutral");
1176
1177 text.setCoordinates(0.85, 0.62, 0.0, Text::CoordinateSystem::NORMAL);
1178
1179 plot.add(text);
1180
1181 text.setColor("00ffff");
1182
1183 text.setText("cation");
1184
1185 text.setCoordinates(0.85, 0.55, 0.0, Text::CoordinateSystem::NORMAL);
1186
1187 plot.add(text);
1188
1189 curve.setXAndY(grid, large);
1190
1191 curve.setColor("00ff00");
1192
1193 curve.setLineWidth(1);
1194
1195 plot.add(curve);
1196
1197 curve.setXAndY(grid, small);
1198
1199 curve.setColor("0000ff");
1200
1201 plot.add(curve);
1202
1203 curve.setXAndY(grid, anion);
1204
1205 curve.setColor("fff00f");
1206
1207 plot.add(curve);
1208
1209 curve.setXAndY(grid, neutral);
1210
1211 curve.setColor("ff00ff");
1212
1213 plot.add(curve);
1214
1215 curve.setXAndY(grid, cation);
1216
1217 curve.setColor("00ffff");
1218
1219 plot.add(curve);
1220
1221 plot.getXAxis().at(0).setTitle(xtitle);
1222
1223 plot.getXAxis().emplace_back();
1224
1225 // plot.getXAxis().at(1).setReciprocalTickFinder();
1226
1227 plot.getXAxis().at(1).setReciprocalLabelFormatter();
1228
1229 plot.getXAxis().at(0).setDrawUnconventionalAxis(false);
1230
1231 plot.getXAxis().at(1).setDrawConventionalAxis(false);
1232
1233 plot.getXAxis().at(1).setDrawConventionalLabels(false);
1234
1235 plot.getXAxis().at(1).setDrawUnconventionalAxis();
1236
1237 plot.getXAxis().at(1).setDrawUnconventionalLabels();
1238
1239 plot.getXAxis().at(1).setTitle(xtitle2);
1240
1241 plot.getYAxis().at(0).setTitle(ytitle);
1242
1243 plots.push_back(plot);
1244
1245 canvas.add(plots);
1246
1247 ostr.str("");
1248
1249 ostr << "Euclidean distance norm: " << spectralfitter.getNorm();
1250
1251 if (std::get<2>(inputdata).size() == std::get<1>(inputdata).size()) {
1252 ostr << " (NNLC; used provided uncertainties)";
1253 } else {
1254 ostr << " (NNLS)";
1255 }
1256
1257 text.setColor("000000");
1258
1259 text.setText(ostr.str());
1260
1261 text.setCoordinates(0.25, 0.8, 0.0);
1262
1263 canvas.add(text);
1264
1265 break;
1266 case Parameters::Arg::CoAdd:
1267 default:
1268
1269 sum = std::vector<double>(spectra.at(0).size());
1270
1271 for (const auto &spectrum : spectra) {
1272
1273 int j = 0;
1274
1275 for (auto s : spectrum) {
1276
1277 sum.at(j++) += s * parameters.getWeights().at(i);
1278 }
1279
1280 ++i;
1281 }
1282
1283 curve.setXAndY(grid, sum);
1284
1285 plot.add(curve);
1286
1287 plot.setYLimits(MinMax::min_max(sum, MinMax::Nice));
1288
1289 plot.getXAxis().at(0).setTitle(xtitle);
1290
1291 plot.getXAxis().emplace_back();
1292
1293 // plot.getXAxis().at(1).setReciprocalTickFinder();
1294
1295 plot.getXAxis().at(1).setReciprocalLabelFormatter();
1296
1297 plot.getXAxis().at(0).setDrawUnconventionalAxis(false);
1298
1299 plot.getXAxis().at(1).setDrawConventionalAxis(false);
1300
1301 plot.getXAxis().at(1).setDrawConventionalLabels(false);
1302
1303 plot.getXAxis().at(1).setDrawUnconventionalAxis();
1304
1305 plot.getXAxis().at(1).setDrawUnconventionalLabels();
1306
1307 plot.getXAxis().at(1).setTitle(xtitle2);
1308
1309 plot.getYAxis().at(0).setTitle(ytitle);
1310
1311 plots.push_back(plot);
1312
1313 canvas.add(plots);
1314
1315 break;
1316 };
1317
1318 std::string f(parameters.getOutputFilename());
1319
1320 FileOutput fileoutput(f.append(".dat"));
1321
1322 spectra.insert(spectra.begin(), grid);
1323
1324 time_t t = time(nullptr);
1325
1326 std::vector<std::string> header{"NASA Ames PAH IR Spectroscopic Database",
1327 std::string("Date: ") + ctime(&t),
1328 "Version : " PROGRAM_VERSION,
1329 "Build: " PROGRAM_BUILD,
1330 std::string("First column: ") + xtitle,
1331 std::string("Other columns: ") + ytitle};
1332
1333 if (parameters.getTool() == Parameters::Arg::SpectralFit) {
1334
1335 header.insert(header.end(), {"First row: UIDs", "Second row: weights"});
1336
1337 std::ostringstream uids;
1338
1339 uids << std::right << std::setw(12) << 0 << '\t';
1340
1341 std::vector<sql_properties> prop = pahdb.getPropertiesByIDs(
1342 parameters.getIds()); // perhaps other name then 'prop'?
1343
1344 for (const auto &properties : prop) {
1345
1346 uids << std::right << std::setw(12) << properties.uid << '\t';
1347 }
1348
1349 header.insert(header.end(), {uids.str(), "First row: weights"});
1350
1351 spectra.at(0).insert(spectra.at(0).begin(), 0);
1352
1353 for (size_t i = 1; i < spectra.size(); i++) {
1354
1355 spectra.at(i).insert(spectra.at(i).begin(),
1356 spectralfitter.getWeights().at(i - 1));
1357 }
1358 }
1359
1360 fileoutput.writeTableHeaderToDisk(header);
1361
1362 fileoutput.writeTableToDisk(spectra);
1363
1364 if (parameters.isX11()) {
1365
1366 canvas.paintOnScreen();
1367 }
1368
1369 if (parameters.isPNG()) {
1370
1371 canvas.paintOnPNG(parameters.getOutputFilename());
1372 }
1373
1374 if (parameters.isJPEG()) {
1375
1376 canvas.paintOnJPEG(parameters.getOutputFilename());
1377 }
1378
1379 if (parameters.isPostscript()) {
1380 canvas.paintOnPostscript(parameters.getOutputFilename());
1381 }
1382
1383 try {
1384
1385 pahdb.setProgress(100.0);
1386 } catch (const Exception &e) {
1387
1388 std::cerr << e.what() << std::endl;
1389
1390 return (1);
1391 }
1392
1393 return (0);
1394}
Curve
Definition Curve.h:11
Line
Definition Line.h:9
PAHdb
Definition PAHdb.h:29
Plot
Definition Plot.h:20
Point
Definition Point.h:11
Text
Definition Text.h:13
From FY2025 onward the NASA Ames PAH IR Spectroscopic Database is being supported through the Laboratory Astrophysics Round 3 directed Work Package at NASA Ames.
From FY2023-2025 the NASA Ames PAH IR Spectroscopic Database was supported through the Laboratory Astrophysics Round 2 directed Work Package at NASA Ames.
From FY2019-2022 the NASA Ames PAH IR Spectroscopic Database was supported through a directed Work Package at NASA Ames titled: "Laboratory Astrophysics – The NASA Ames PAH IR Spectroscopic Database".
© Copyright 2021-2025, Christiaan Boersma