From 95cd9ddad32fba1f0b8a94311c62e1ee12d27cc1 Mon Sep 17 00:00:00 2001 From: Steven Haemmerich Date: Mon, 27 Sep 2021 15:27:16 +0200 Subject: [PATCH 1/2] stuff for require gsl --- src/misc/astronomy/hardnessratio.sl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/misc/astronomy/hardnessratio.sl b/src/misc/astronomy/hardnessratio.sl index d0d085b9..dff4fc52 100644 --- a/src/misc/astronomy/hardnessratio.sl +++ b/src/misc/astronomy/hardnessratio.sl @@ -1,3 +1,5 @@ +require( "gsl" ); + define behr() %!%+ %\function{behr} @@ -5,7 +7,7 @@ define behr() %\usage{Double_Type (HR , HR_err_max , HR_err_min) = behr(Double_Type hard_count, Double_Type soft_count);} %\qualifiers{ % \qualifier{scale_s}{: Scaling factor for soft_count; soft_count = scale_s*soft_count (default = 1)} -% \qualifier{scale_h}{: Scaling factor for hard_count; hard_count = hard_s*soft_count (default = 1)} +% \qualifier{scale_h}{: Scaling factor for hard_count; hard_count = scale_h*hard_count (default = 1)} % \qualifier{back_s}{: background counts in the soft band} % \qualifier{back_h}{: background counts in the hard band} % \qualifier{bkg_scale_s}{: Scaling factor for source background to observed background counts in the soft band. @@ -30,7 +32,7 @@ define behr() % Use qualifier if also the grid used (2000 color values ranging from -0.999 to 0.999), and the % normalized probability distribution are desired. % -% Requires gsl, make sure module is loaded! Based on Park et al., 2006, ApJ, 652, 610. Based on original isis-code by Mike Nowak. +% Based on Park et al., 2006, ApJ, 652, 610. Based on original isis-code by Mike Nowak. % %\seealso{hardnessratio;} %!%- @@ -101,14 +103,14 @@ define behr() variable iwa = where(phsum <= pval/2.); variable iwb = where(phsum >= 1-pval/2.); if ((length(bhr[iwb]) == 0) or (length(bhr[iwa]) == 0)){ - print(sprintf("%s: Ranges for uncertainties are empty. Make sure you have run require('gsl') before!",_function_name())); + print(sprintf("%s: Ranges for uncertainties are empty.",_function_name())); print(sprintf("%s: Will only return best value and _NaN for uncertanties.",_function_name())); return (bhr_max[0] , _NaN , _NaN); } variable bhrmin = max(bhr[iwa]); variable bhrmax = min(bhr[iwb]); if ((bhr_max[0] <= bhrmin[0]) or (bhr_max[0] >= bhrmax[0])){ - print(sprintf("%s: Uncertainty-range does not agree with best value. Make sure you have run require('gsl') before!",_function_name())); + print(sprintf("%s: Uncertainty-range does not agree with best value.",_function_name())); print(sprintf("%s: Will only return best value and _NaN for uncertanties.",_function_name())); return (bhr_max[0] , _NaN , _NaN); } -- GitLab From c641d7be0ef8ca22f7fb4f765097d8da77fbe6c3 Mon Sep 17 00:00:00 2001 From: Steven Haemmerich Date: Mon, 27 Sep 2021 16:45:54 +0200 Subject: [PATCH 2/2] fix behr function; now gsl-stuff actually works --- src/misc/astronomy/hardnessratio.sl | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/misc/astronomy/hardnessratio.sl b/src/misc/astronomy/hardnessratio.sl index dff4fc52..c4c1d81e 100644 --- a/src/misc/astronomy/hardnessratio.sl +++ b/src/misc/astronomy/hardnessratio.sl @@ -1,5 +1,8 @@ -require( "gsl" ); - +require("gsl" , "gsl"); +% don't ask me why but it wont work otherwise. The machine knows! - M. Scott +% +% +% define behr() %!%+ %\function{behr} @@ -80,9 +83,9 @@ define behr() max_k = min([int(b) , n - j]); } _for k (0,max_k,1){ - pjk = lngamma(da + j + 1) + lngamma(db + k + 1) - lngamma(da + 1) - lngamma(db + 1) - - j*log(ga + 1.) - k*log(gb + 1) + lngamma(a + b + 2 - j - k) - - lngamma(a - j + 1) - lngamma(b - k + 1) - lngamma(j + 1) - lngamma(k + 1) + pjk = gsl->lngamma(da + j + 1) + gsl->lngamma(db + k + 1) - gsl->lngamma(da + 1) - gsl->lngamma(db + 1) + - j*log(ga + 1.) - k*log(gb + 1) + gsl->lngamma(a + b + 2 - j - k) + - gsl->lngamma(a - j + 1) - gsl->lngamma(b - k + 1) - gsl->lngamma(j + 1) - gsl->lngamma(k + 1) + (a - j)*log(1 + bhr) + (b - k)*log(1 - bhr) - (a + b + 2 - j - k)*log(ea + eb + (ea - eb)*bhr); psum += exp(__tmp(pjk)); @@ -91,7 +94,7 @@ define behr() psum = psum/sum(psum*dbhr); } else{ - variable ph = lngamma(a + b + 2) + log(2.) - lngamma(a + 1) - lngamma(b + 1) + variable ph = gsl->lngamma(a + b + 2) + log(2.) - gsl->lngamma(a + 1) - gsl->lngamma(b + 1) + (a + 1)*log(ea) + (b + 1)*log(eb) + a*log(1 + bhr) + b*log(1 - bhr) - (a + b + 2)*log(ea + eb + (ea - eb)*bhr); psum = exp(__tmp(ph)); -- GitLab