Add mass_profile tools

* These tools are mainly use to calculate the total gravitational mass
  profile, as well as the intermediate products (e.g., surface
  brightness profile fitting, gas density profile, NFW fitting, etc.)
* There are additional tools for calculating the luminosity and flux.
* These tools mainly developed by Junhua GU, and contributed by
  Weitian (Aaron) LI, and Zhenghao ZHU.

1 files changed, 124 insertions, 0 deletions

diff --git a/mass_profile/adapt_trapezoid.h b/mass_profile/adapt_trapezoid.h
new file mode 100644
index 0000000..35d22b4
--- /dev/null
+++ b/mass_profile/adapt_trapezoid.h
@@ -0,0 +1,124 @@
+#ifndef ADAPT_TRAPEZOID_H
+#define ADAPT_TRAPEZOID_H
+#include <list>
+#include <utility>
+#include <cassert>
+//#include <iostream>
+
+//using namespace std;
+
+template<typename T1,typename T2,typename T3>
+class triple
+{
+ public:
+  T1 first;
+  T2 second;
+  T3 third;
+  triple(T1 x1,T2 x2,T3 x3)
+    :first(x1),second(x2),third(x3)
+    {
+    }
+};
+
+
+template<typename T1,typename T2,typename T3>
+triple<T1,T2,T3> make_triple(T1 x1,T2 x2,T3 x3)
+{
+  return triple<T1,T2,T3>(x1,x2,x3);
+}
+
+template <typename T>
+T trapezoid(T (*fun)(T),T x1,T x2,T err_limit)
+{
+  int n=256;
+  T result;
+  const int max_division=24;
+  T old_value=0;
+  for(int i=1;i<max_division;i++)
+    {
+      result=0.;
+      n*=2;
+      T step=(x2-x1)/n;
+      for(int j=0;j<n;j++)
+	{
+	  result+=(fun(x1+(j+1)*step)+fun(x1+j*step))*step/T(2.);
+	}
+      old_value-=result;
+      old_value=old_value<0?-old_value:old_value;
+      if(old_value<err_limit)
+	{
+	  return result;
+	}
+      old_value=result;
+    }
+}
+
+
+template <typename T>
+T adapt_trapezoid(T (*fun)(T),T x1,T x2,T err_limit)
+{
+  //  const err_limit=.001;
+  typedef triple<T,T,bool> interval;
+  /*T for interval type,
+    bool for state trur for still to be updated,
+    false for do not need to be updated
+   */
+  std::list<interval> interval_list;
+  T current_sum=((fun(x1)+fun(x2))/2.*(x2-x1));
+  interval_list.push_back(make_triple(x1,current_sum,true));
+  interval_list.push_back(make_triple(x2,(T)0.,true));
+  bool int_state=1;
+  int n_intervals=1;
+  while(int_state)
+    {
+      //std::cout<<n_intervals<<std::endl;
+      int_state=0;
+      typename std::list<interval>::iterator i1=interval_list.begin();
+      typename std::list<interval>::iterator i2=interval_list.begin();
+      i2++;
+      for(;i2!=interval_list.end();i1++,i2=i1,i2++)
+	{
+	  //cout<<i1->first<<"\t"<<i2->first<<endl;
+	  //assert(i2->first>i1->first);
+	  if(i1->third)
+	    {
+	      interval new_interval((i1->first+i2->first)/2,0,true);
+	      
+	      T sum1,sum2;
+	      sum1=(fun(new_interval.first)+fun(i1->first))/2*(new_interval.first-i1->first);
+	      sum2=(fun(new_interval.first)+fun(i2->first))/2*(i2->first-new_interval.first);
+	      new_interval.second=sum2;
+	      T err;
+	      err=i1->second-sum1-sum2;
+	      err=err<0?-err:err;
+
+	      if(err>err_limit/n_intervals)
+		{
+		  i1->second=sum1;
+		  interval_list.insert(i2,new_interval);
+		  n_intervals++;
+		  if(n_intervals>10e6)
+		    {
+		      
+		      break;
+		    }
+		}
+	      else
+		{
+		  i1->third=false;
+		}
+	      int_state=1;
+	    }
+	}
+      
+    }
+  T result=0;
+  for(typename std::list<interval>::iterator i=interval_list.begin();i!=interval_list.end();i++)
+    {
+      result+=i->second;
+    }
+  return result;
+}
+
+#endif
+//end of the file