double Funktionswertmu(int mu,double x,double y,double z,double** &exponenten,double** &Coordinates,int** &L,double** &contraccoeff,int* &mubasisfunktionentyp,int* &mucoordinaten,int* &prim){
	double WERT=0.0;
	for(int p=0; p<prim[mubasisfunktionentyp[mu]]; p++){
		WERT=WERT+contraccoeff[mubasisfunktionentyp[mu]][p]*pow((x-Coordinates[mucoordinaten[mu]][0]),L[mubasisfunktionentyp[mu]][0])*
				pow((y-Coordinates[mucoordinaten[mu]][1]),L[mubasisfunktionentyp[mu]][1])*
				pow((z-Coordinates[mucoordinaten[mu]][2]),L[mubasisfunktionentyp[mu]][2])*
				exp(-exponenten[mubasisfunktionentyp[mu]][p]*pow((x-Coordinates[mucoordinaten[mu]][0]),2))*
				exp(-exponenten[mubasisfunktionentyp[mu]][p]*pow((y-Coordinates[mucoordinaten[mu]][1]),2))*
				exp(-exponenten[mubasisfunktionentyp[mu]][p]*pow((z-Coordinates[mucoordinaten[mu]][2]),2));
	}
	return WERT;
}

double CalcRHOXYZ(double x,double y,double z,double** &exponenten,double** &Coordinates,int** &L,double** &contraccoeff,int* &mubasisfunktionentyp,int* &mucoordinaten,double* &Pdiff,int* &prim,int N){
double RHO=0.0;
	for(int mu=0; mu<N; mu++){
		for(int nu=0; nu<N; nu++){
			RHO+=Pdiff[mu+nu*N]*Funktionswertmu(mu,x,y,z,exponenten,Coordinates,L,contraccoeff,mubasisfunktionentyp,mucoordinaten,prim)*
					Funktionswertmu(nu,x,y,z,exponenten,Coordinates,L,contraccoeff,mubasisfunktionentyp,mucoordinaten,prim);
		}
	}
 return RHO;
}


// Berechnung der Elektronendichte RHO am Ort x,y,z.

	//test
    //cout<<CalcRHOXYZ(1.0,1.0,1.0,exponenten,Coordinates,L,contraccoeff,mubasisfunktionentyp,mucoordinaten,Pdiff,prim,N);

	// Berechnung des RHOarrays

	// Die Grenzen eingeben,zwischen denen abgerastert werden soll
	double XMax=-1.0;
	double XMin=-3.0;
	double YMax=4.0;
	double YMin=2.0;
	double ZMax=1.0;
	double ZMin=-1.0;
	// Die Anzahl der Rasterpunkte eingeben in Richtung X,Y,Z
	int MaxPointsX=20;
	int MaxPointsY=25;
	int MaxPointsZ=30;
	// Berechnung der Schrittweite zwischen den Rasterpunkten
	double SchrittweiteX=abs(XMin-XMax)/(MaxPointsX-1);
	double SchrittweiteY=abs(YMin-YMax)/(MaxPointsY-1);
	double SchrittweiteZ=abs(ZMin-ZMax)/(MaxPointsZ-1);
	// Festlegung des Arrays in dem die Elektronendichten zu einem Rasterpunkt eingeschrieben werden
	double* RHOARRAY = new double [MaxPointsX*MaxPointsY*MaxPointsZ];
	// Festlegung, welche Achse konstant gehalten wird. z.B. für E-dichte in XY-Ebene z-Achse konstant halten
	bool constantXAchse=true;
	bool constantYAchse=false;
	bool constantZAchse=false;
	// Festlegung wo auf der konstanten Achse der Querschnitt erfolgt
	double Xconstant=1.345;
	double Yconstant=1.2392;
	double Zconstant=1.782346;

	int zaehler=0;
	if(constantXAchse==true){MaxPointsX=1;}
	if(constantYAchse==true){MaxPointsY=1;}
	if(constantZAchse==true){MaxPointsZ=1;}
	for(int PointX=0; PointX<MaxPointsX; PointX++){
		for(int PointY=0; PointY<MaxPointsY; PointY++){
			for(int PointZ=0; PointZ<MaxPointsZ; PointZ++){
				double x=XMin+SchrittweiteX*PointX;
				double y=YMin+SchrittweiteY*PointY;
				double z=ZMin+SchrittweiteZ*PointZ;
				if(constantXAchse==true){x=Xconstant;}
				if(constantYAchse==true){y=Yconstant;}
				if(constantZAchse==true){z=Zconstant;}
				RHOARRAY[PointX+PointY*MaxPointsX+PointZ*MaxPointsX*MaxPointsY]=CalcRHOXYZ(x,y,z,exponenten,Coordinates,L,contraccoeff,mubasisfunktionentyp,mucoordinaten,Pdiff,prim,N);
				zaehler+=1;
			}
		}
	}


	// Ausgabe der Elektronendichte an den Rasterpunkten und deren Koordinaten
	cout<<endl;
	cout<<"Ausgabe der Elektronendichte an den Rasterpunkten"<<endl;
	cout<<endl;

	cout<<endl;
	cout<<"Anzahl berechneter Rasterpunkte für die Elektronendichte-Darstellung"<<endl;
	cout<<zaehler<<endl;
	cout<<endl;

	if(constantXAchse==false){cout<<" x-Koordinaten          ";}
	if(constantYAchse==false){cout<<"y-Kooridnaten          ";}
	if(constantZAchse==false){cout<<"z-Koordinaten           ";}
	cout<<"   DELTARHO(x,y,z)      ";
	if(constantXAchse==true){printf("x-Querschn.:%12.8f" ,Xconstant);}
	if(constantYAchse==true){printf("  y-Querschn.:%12.8f" ,Yconstant);}
	if(constantZAchse==true){printf("  z-Querschn.:%12.8f" ,Zconstant);}
	cout<<endl;
	cout<<endl;

	for(int PointX=0; PointX<MaxPointsX; PointX++){
		for(int PointY=0; PointY<MaxPointsY; PointY++){
			for(int PointZ=0; PointZ<MaxPointsZ; PointZ++){
				double x=XMin+SchrittweiteX*PointX;
				double y=YMin+SchrittweiteY*PointY;
				double z=ZMin+SchrittweiteZ*PointZ;
				if(constantXAchse==false){printf("%12.8f           ", x);}
				if(constantYAchse==false){printf("%12.8f           ", y);}
				if(constantZAchse==false){printf("%12.8f           ", z);}
				cout<<"      ";
				cout<<RHOARRAY[PointX+PointY*MaxPointsX+PointZ*MaxPointsX*MaxPointsY]<<endl;
			}
		}
	}