From 34ed04686523ecdafb96c6cdff379df176cf1984 Mon Sep 17 00:00:00 2001 From: Aaron LI Date: Sat, 28 May 2016 09:40:25 +0800 Subject: Add astro/cosmo_calc from "chandra-acis-analysis" --- astro/cosmo_calc/ddivid.cc | 56 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 56 insertions(+) create mode 100644 astro/cosmo_calc/ddivid.cc (limited to 'astro/cosmo_calc/ddivid.cc') diff --git a/astro/cosmo_calc/ddivid.cc b/astro/cosmo_calc/ddivid.cc new file mode 100644 index 0000000..aa72471 --- /dev/null +++ b/astro/cosmo_calc/ddivid.cc @@ -0,0 +1,56 @@ +#include + +double ddivid(double (*foo)(double),double x1,double x2,double err) +{ + assert(x2>x1); + assert(foo(x2)*foo(x1)<0); + while(x2-x1>err) + { + double x=(x2+x1)/2.; + double y=foo(x); + double y1=foo(x1); + double y2=foo(x2); + if(y1*y<0) + { + x2=x; + } + else if(y2*y<0) + { + x1=x; + } + else + { + assert(0); + } + + } + return (x1+x2)/2.; +} + + +double ddivid(double (*foo)(double),double z,double x1,double x2,double err) +{ + assert(x2>x1); + assert((foo(x2)-z)*(foo(x1)-z)<0); + while(x2-x1>err) + { + double x=(x2+x1)/2.; + double y=foo(x)-z; + double y1=foo(x1)-z; + double y2=foo(x2)-z; + if(y1*y<0) + { + x2=x; + } + else if(y2*y<0) + { + x1=x; + } + else + { + assert(0); + } + + } + return (x1+x2)/2.; +} -- cgit v1.2.2