////////// Trapezregel (sechs Parameter) 
 
 
#define FUNC(z, a, b, rho12, rho13, rho23) ((*func)(z, a, b, rho12, rho13, rho23))
 
double trapzd6(double (*func)(double z, double, double, double, double, double),double a, double b, double rho12, double rho13, double rho23, double unteregrenze, double oberegrenze, int n)
 
{
double sum,h;
static double s;
//static double s, a, b, rho12, rho13, rho23;
int it;
if (n == 1) 
	{
	return s = 0.5 * (oberegrenze - unteregrenze) * 
		(FUNC(unteregrenze, a, b, rho12, rho13, rho23) 
		+ FUNC(oberegrenze, a, b, rho12, rho13, rho23));
	} 
else 
	{ h = (oberegrenze - unteregrenze)/n;
	  sum = 0.5 * (FUNC(unteregrenze, a, b, rho12,  rho13,rho23)    + FUNC(oberegrenze, a, b, rho12, rho13, rho23));
	  for (it=1;it<n;it++) 
	  {sum = sum + FUNC(unteregrenze + it*h, a, b, rho12, rho13, rho23);
	  }
 
	return h*sum;
	}
}
 
 
////////// Hilfsfunktion für N3 (für Integralreduktion) 
 
 
double hilfsfunktion(double z, double a, double b, double rho12, double rho13, double rho23)
{
return cnd(z) * cbnd((b-rho23*z)/(sqrt(1-rho23 * rho23)),(a-rho13*z)/(sqrt(1- rho13 * rho13)), (rho12-rho13*rho23)/((sqrt(1- rho13 * rho13))*(sqrt(1- rho23 * rho23))));
 
}
// cnd und cbnd hab ich vorher programmiert, die klappen.
 
////////// ctnd-Funktion
 
double ctnd(double a, double b, double c, double rho12, double rho13, double rho23)
{
 
return trapzd6(hilfsfunktion,a,b,rho12, rho13, rho23, -10/*untergrenze*/, c/*oberegrenze*/, 100000 /*Schritte*/);
 
}
 
 
////////////////////////// main
 
void main (void)
{ 
 
 
printf("\n\n Der Wert ist %lf", ctnd(0, 0, 0, 0.5, 0.5, 0.5));
 
 
}
 

////////// Hilfsfunktion g (für Integralreduktion) 
 
double g(double z, double a, double b)
{
return exp(a*z)+exp(b*z);
}
 
//g hab ich jetzt einfach mal so genannt, in Wirklichkeit ist sie komplizierter.
 
////////// Funktion f
 
double f(double a, double b, double c)
{
return trapzd(g, -20.0, c, 10.0);
 
}
 
////////////////////////// main
 
void main (void)
{ 
 
 
printf("\n\n Der Wert ist %lf", f(0.3, 0.4, 0.5));
}

#define FUNC(x) ((*func)(x))
double trapzd(double (*func)(double), double a, double b, int n)
//This routine computes the nth stage of refinement of an extended trapezoidal rule. func is input
//as a pointer to the function to be integrated between limits a and b, also input. When called with
//n=1, the routine returns the crudest estimate of .b
//a f(x)dx. Subsequent calls with n=2,3,...
//(in that sequential order) will improve the accuracy by adding 2n-2 additional interior points.
{
double x,tnm,sum,del;
static double s;
int it,j;
if (n == 1) {
return (s=0.5*(b-a)*(FUNC(a)+FUNC(b)));
} else {
for (it=1,j=1;j<n-1;j++) it <<= 1;
tnm=it;
del=(b-a)/tnm; //This is the spacing of the points to be added.
x=a+0.5*del;
for (sum=0.0,j=1;j<=it;j++,x+=del) sum += FUNC(x);
s=0.5*(s+(b-a)*sum/tnm); //This replaces s by its refined value.
return s;
}
}