/* turbulence model: RSM wall reflection model by Craft and Launder (1992), (see Paper Craft T.J., Launder B.E., 2001)*/
/* UDF (User Defined Function) Routine for FLUENT                                                             */
/* see also the paper:                                                                                        */
/* Implementation and comparison of different turbulence models for three dimensional wall jets with Fluent   */
/* FLUENT CFD Forum 2005 Bad Nauheim, Germany								      */
/* Heschl Ch*., Sanz W.**, Klanatsky P.*                                                                      */
/* *  Fachhochschulstudiengänge Burgenland GmbH, Steinamangerstraße 21, A-7423 Pinkafeld                      */
/* ** TU Graz, Institut of Thermal Turbomachinery and Machine Dynamics, Inffeldgasse 25, A-8010 Graz          */
/* Contact: christian.heschl@fh-burgenland.at, wolfgang.sanz@tugraz.at, peter.klanatsky@fh-burgenland.at      */

/* Die Autoren übernhemen keine Verantwortung oder Haftung für den Inhalt der unten angeführten UDF-Routinen  */



#include "udf.h"
#include "math.h"

/* Turbulence model constants */
const real c1s=0.5;
const real c2s=0.08;
const real c2as=0.1;
const real c2bs=0.4;
const real cl=2.55;



DEFINE_SOURCE(uu_source, c, t, dS, eqn)
{
real source;
real u11, u12, u13, u21, u22, u23, u31, u32, u33;
real a11, a12, a13, a21, a22, a23, a31, a32, a33;
real n1, n2, n3, f;
real phi_11_w1, phi_11_w2;

f=pow(C_K(c,t),1.5)/C_D(c,t)/cl/C_UDSI(c,t,0);


n1=C_UDMI(c,t,0);
n2=C_UDMI(c,t,1);
n3=C_UDMI(c,t,2);


u11=C_RUU(c,t);
u12=C_RUV(c,t);
u13=C_RUW(c,t);
u21=C_RUV(c,t);
u22=C_RVV(c,t);
u23=C_RVW(c,t);
u31=C_RUW(c,t);
u32=C_RVW(c,t);
u33=C_RWW(c,t);

a11=u11/C_K(c,t)-2./3.;
a12=u12/C_K(c,t);
a13=u13/C_K(c,t);
a21=u21/C_K(c,t);
a22=u22/C_K(c,t)-2./3.;
a23=u23/C_K(c,t);
a31=u31/C_K(c,t);
a32=u32/C_K(c,t);
a33=u33/C_K(c,t)-2./3.;



phi_11_w1=c1s*C_D(c,t)/C_K(c,t)*(
								 u11*n1*n1+u12*n1*n2+u13*n1*n3+u21*n2*n1+u22*n2*n2+u23*n2*n3+u31*n3*n1+u32*n3*n2+u33*n3*n3
								 -3./2.*(u11*n1*n1+u21*n2*n1+u31*n3*n1)
								 -3./2.*(u11*n1*n1+u21*n1*n2+u31*n1*n3)
								 )*f;
phi_11_w2=-c2s*(
				C_DUDX(c,t)*u11+C_DVDX(c,t)*u21+C_DWDX(c,t)*u31+
				C_DUDY(c,t)*u12+C_DVDY(c,t)*u22+C_DWDY(c,t)*u32+
				C_DUDZ(c,t)*u13+C_DVDZ(c,t)*u23+C_DWDZ(c,t)*u33
				)*(n1*n1+n2*n2+n3*n3-3.*n1*n1)*f;

phi_11_w2=phi_11_w2-
			c2as*C_K(c,t)*(

			C_DUDX(c,t)*n1*n1*a11+C_DUDY(c,t)*n1*n1*a12+C_DUDZ(c,t)*n1*n1*a13+
			C_DVDX(c,t)*n1*n2*a11+C_DVDY(c,t)*n1*n2*a12+C_DVDZ(c,t)*n1*n2*a13+
			C_DWDX(c,t)*n1*n3*a11+C_DWDY(c,t)*n1*n3*a12+C_DWDZ(c,t)*n1*n3*a13+

			C_DUDX(c,t)*n2*n1*a21+C_DUDY(c,t)*n2*n1*a22+C_DUDZ(c,t)*n2*n1*a23+
			C_DVDX(c,t)*n2*n2*a21+C_DVDY(c,t)*n2*n2*a22+C_DVDZ(c,t)*n2*n2*a23+
			C_DWDX(c,t)*n2*n3*a21+C_DWDY(c,t)*n2*n3*a22+C_DWDZ(c,t)*n2*n3*a23+

			C_DUDX(c,t)*n3*n1*a31+C_DUDY(c,t)*n3*n1*a32+C_DUDZ(c,t)*n3*n1*a33+
			C_DVDX(c,t)*n3*n2*a31+C_DVDY(c,t)*n3*n2*a32+C_DVDZ(c,t)*n3*n2*a33+
			C_DWDX(c,t)*n3*n3*a31+C_DWDY(c,t)*n3*n3*a32+C_DWDZ(c,t)*n3*n3*a33

			-3./2.*(C_DUDX(c,t)*n1*n1*a11+C_DUDY(c,t)*n1*n1*a12+C_DUDZ(c,t)*n1*n1*a13+
					C_DUDX(c,t)*n2*n1*a21+C_DUDY(c,t)*n2*n1*a22+C_DUDZ(c,t)*n2*n1*a23+
					C_DUDX(c,t)*n3*n1*a31+C_DUDY(c,t)*n3*n1*a32+C_DUDZ(c,t)*n3*n1*a33)

			-3./2.*(C_DUDX(c,t)*n1*n1*a11+C_DUDY(c,t)*n1*n1*a12+C_DUDZ(c,t)*n1*n1*a13+
					C_DUDX(c,t)*n2*n1*a21+C_DUDY(c,t)*n2*n1*a22+C_DUDZ(c,t)*n2*n1*a23+
					C_DUDX(c,t)*n3*n1*a31+C_DUDY(c,t)*n3*n1*a32+C_DUDZ(c,t)*n3*n1*a33)
							
						)*f;


phi_11_w2=phi_11_w2+
			c2bs*C_K(c,t)*(
			
			C_DUDX(c,t)*n1*n1+C_DUDY(c,t)*n1*n2+C_DUDZ(c,t)*n1*n3+
			C_DVDX(c,t)*n2*n1+C_DVDY(c,t)*n2*n2+C_DVDZ(c,t)*n2*n3+
			C_DWDX(c,t)*n3*n1+C_DWDY(c,t)*n3*n2+C_DWDZ(c,t)*n3*n3
		
						)

			*(n1*n1-1./3.*(n1*n1+n2*n2+n3*n3))*f;


// prüfen ob die Dichte berücksichtigt werden muss!!

	dS[eqn]= -c1s*C_D(c,t)/C_K(c,t)*2.*n1*n1*f;
	source= phi_11_w1+phi_11_w2;

return source;
}


DEFINE_SOURCE(vv_source, c, t, dS, eqn)
{
real source;
real u11, u12, u13, u21, u22, u23, u31, u32, u33;
real a11, a12, a13, a21, a22, a23, a31, a32, a33;
real n1, n2, n3, f;
real phi_22_w1, phi_22_w2;

f=pow(C_K(c,t),1.5)/C_D(c,t)/cl/C_UDSI(c,t,0);


n1=C_UDMI(c,t,0);
n2=C_UDMI(c,t,1);
n3=C_UDMI(c,t,2);


u11=C_RUU(c,t);
u12=C_RUV(c,t);
u13=C_RUW(c,t);
u21=C_RUV(c,t);
u22=C_RVV(c,t);
u23=C_RVW(c,t);
u31=C_RUW(c,t);
u32=C_RVW(c,t);
u33=C_RWW(c,t);

a11=u11/C_K(c,t)-2./3.;
a12=u12/C_K(c,t);
a13=u13/C_K(c,t);
a21=u21/C_K(c,t);
a22=u22/C_K(c,t)-2./3.;
a23=u23/C_K(c,t);
a31=u31/C_K(c,t);
a32=u32/C_K(c,t);
a33=u33/C_K(c,t)-2./3.;



phi_22_w1=c1s*C_D(c,t)/C_K(c,t)*(
								 u11*n1*n1+u12*n1*n2+u13*n1*n3+u21*n2*n1+u22*n2*n2+u23*n2*n3+u31*n3*n1+u32*n3*n2+u33*n3*n3
								 -3./2.*(u12*n1*n2+u22*n2*n2+u32*n3*n2)
								 -3./2.*(u12*n2*n1+u22*n2*n2+u32*n2*n3)
								 )*f;
phi_22_w2=-c2s*(
				C_DUDX(c,t)*u11+C_DVDX(c,t)*u21+C_DWDX(c,t)*u31+
				C_DUDY(c,t)*u12+C_DVDY(c,t)*u22+C_DWDY(c,t)*u32+
				C_DUDZ(c,t)*u13+C_DVDZ(c,t)*u23+C_DWDZ(c,t)*u33
				)*(n1*n1+n2*n2+n3*n3-3.*n2*n2)*f;

phi_22_w2=phi_22_w2-
			c2as*C_K(c,t)*(

			C_DUDX(c,t)*n1*n1*a11+C_DUDY(c,t)*n1*n1*a12+C_DUDZ(c,t)*n1*n1*a13+
			C_DVDX(c,t)*n1*n2*a11+C_DVDY(c,t)*n1*n2*a12+C_DVDZ(c,t)*n1*n2*a13+
			C_DWDX(c,t)*n1*n3*a11+C_DWDY(c,t)*n1*n3*a12+C_DWDZ(c,t)*n1*n3*a13+

			C_DUDX(c,t)*n2*n1*a21+C_DUDY(c,t)*n2*n1*a22+C_DUDZ(c,t)*n2*n1*a23+
			C_DVDX(c,t)*n2*n2*a21+C_DVDY(c,t)*n2*n2*a22+C_DVDZ(c,t)*n2*n2*a23+
			C_DWDX(c,t)*n2*n3*a21+C_DWDY(c,t)*n2*n3*a22+C_DWDZ(c,t)*n2*n3*a23+

			C_DUDX(c,t)*n3*n1*a31+C_DUDY(c,t)*n3*n1*a32+C_DUDZ(c,t)*n3*n1*a33+
			C_DVDX(c,t)*n3*n2*a31+C_DVDY(c,t)*n3*n2*a32+C_DVDZ(c,t)*n3*n2*a33+
			C_DWDX(c,t)*n3*n3*a31+C_DWDY(c,t)*n3*n3*a32+C_DWDZ(c,t)*n3*n3*a33

			-3./2.*(C_DVDX(c,t)*n1*n2*a11+C_DVDY(c,t)*n1*n2*a12+C_DVDZ(c,t)*n1*n2*a13+
					C_DVDX(c,t)*n2*n2*a21+C_DVDY(c,t)*n2*n2*a22+C_DVDZ(c,t)*n2*n2*a23+
					C_DVDX(c,t)*n3*n2*a31+C_DVDY(c,t)*n3*n2*a32+C_DVDZ(c,t)*n3*n2*a33)

			-3./2.*(C_DVDX(c,t)*n1*n2*a11+C_DVDY(c,t)*n1*n2*a12+C_DVDZ(c,t)*n1*n2*a13+
					C_DVDX(c,t)*n2*n2*a21+C_DVDY(c,t)*n2*n2*a22+C_DVDZ(c,t)*n2*n2*a23+
					C_DVDX(c,t)*n3*n2*a31+C_DVDY(c,t)*n3*n2*a32+C_DVDZ(c,t)*n3*n2*a33)
							
						)*f;


phi_22_w2=phi_22_w2+
			c2bs*C_K(c,t)*(
			
			C_DUDX(c,t)*n1*n1+C_DUDY(c,t)*n1*n2+C_DUDZ(c,t)*n1*n3+
			C_DVDX(c,t)*n2*n1+C_DVDY(c,t)*n2*n2+C_DVDZ(c,t)*n2*n3+
			C_DWDX(c,t)*n3*n1+C_DWDY(c,t)*n3*n2+C_DWDZ(c,t)*n3*n3
		
						)

			*(n2*n2-1./3.*(n1*n1+n2*n2+n3*n3))*f;


// prüfen ob die Dichte berücksichtigt werden muss!!

	dS[eqn]= -c1s*C_D(c,t)/C_K(c,t)*2.*n2*n2*f;
	source= phi_22_w1+phi_22_w2;

return source;
}


DEFINE_SOURCE(ww_source, c, t, dS, eqn)
{
real source;
real u11, u12, u13, u21, u22, u23, u31, u32, u33;
real a11, a12, a13, a21, a22, a23, a31, a32, a33;
real n1, n2, n3, f;
real phi_33_w1, phi_33_w2;

f=pow(C_K(c,t),1.5)/C_D(c,t)/cl/C_UDSI(c,t,0);


n1=C_UDMI(c,t,0);
n2=C_UDMI(c,t,1);
n3=C_UDMI(c,t,2);


u11=C_RUU(c,t);
u12=C_RUV(c,t);
u13=C_RUW(c,t);
u21=C_RUV(c,t);
u22=C_RVV(c,t);
u23=C_RVW(c,t);
u31=C_RUW(c,t);
u32=C_RVW(c,t);
u33=C_RWW(c,t);

a11=u11/C_K(c,t)-2./3.;
a12=u12/C_K(c,t);
a13=u13/C_K(c,t);
a21=u21/C_K(c,t);
a22=u22/C_K(c,t)-2./3.;
a23=u23/C_K(c,t);
a31=u31/C_K(c,t);
a32=u32/C_K(c,t);
a33=u33/C_K(c,t)-2./3.;



phi_33_w1=c1s*C_D(c,t)/C_K(c,t)*(
								 u11*n1*n1+u12*n1*n2+u13*n1*n3+u21*n2*n1+u22*n2*n2+u23*n2*n3+u31*n3*n1+u32*n3*n2+u33*n3*n3
								 -3./2.*(u13*n1*n3+u23*n2*n3+u33*n3*n3)
								 -3./2.*(u13*n3*n1+u23*n3*n2+u33*n3*n3)
								 )*f;
phi_33_w2=-c2s*(
				C_DUDX(c,t)*u11+C_DVDX(c,t)*u21+C_DWDX(c,t)*u31+
				C_DUDY(c,t)*u12+C_DVDY(c,t)*u22+C_DWDY(c,t)*u32+
				C_DUDZ(c,t)*u13+C_DVDZ(c,t)*u23+C_DWDZ(c,t)*u33
				)*(n1*n1+n2*n2+n3*n3-3.*n3*n3)*f;

phi_33_w2=phi_33_w2-
			c2as*C_K(c,t)*(

			C_DUDX(c,t)*n1*n1*a11+C_DUDY(c,t)*n1*n1*a12+C_DUDZ(c,t)*n1*n1*a13+
			C_DVDX(c,t)*n1*n2*a11+C_DVDY(c,t)*n1*n2*a12+C_DVDZ(c,t)*n1*n2*a13+
			C_DWDX(c,t)*n1*n3*a11+C_DWDY(c,t)*n1*n3*a12+C_DWDZ(c,t)*n1*n3*a13+

			C_DUDX(c,t)*n2*n1*a21+C_DUDY(c,t)*n2*n1*a22+C_DUDZ(c,t)*n2*n1*a23+
			C_DVDX(c,t)*n2*n2*a21+C_DVDY(c,t)*n2*n2*a22+C_DVDZ(c,t)*n2*n2*a23+
			C_DWDX(c,t)*n2*n3*a21+C_DWDY(c,t)*n2*n3*a22+C_DWDZ(c,t)*n2*n3*a23+

			C_DUDX(c,t)*n3*n1*a31+C_DUDY(c,t)*n3*n1*a32+C_DUDZ(c,t)*n3*n1*a33+
			C_DVDX(c,t)*n3*n2*a31+C_DVDY(c,t)*n3*n2*a32+C_DVDZ(c,t)*n3*n2*a33+
			C_DWDX(c,t)*n3*n3*a31+C_DWDY(c,t)*n3*n3*a32+C_DWDZ(c,t)*n3*n3*a33

			-3./2.*(C_DWDX(c,t)*n1*n3*a11+C_DWDY(c,t)*n1*n3*a12+C_DWDZ(c,t)*n1*n3*a13+
					C_DWDX(c,t)*n2*n3*a21+C_DWDY(c,t)*n2*n3*a22+C_DWDZ(c,t)*n2*n3*a23+
					C_DWDX(c,t)*n3*n3*a31+C_DWDY(c,t)*n3*n3*a32+C_DWDZ(c,t)*n3*n3*a33)

			-3./2.*(C_DWDX(c,t)*n1*n3*a11+C_DWDY(c,t)*n1*n3*a12+C_DWDZ(c,t)*n1*n3*a13+
					C_DWDX(c,t)*n2*n3*a21+C_DWDY(c,t)*n2*n3*a22+C_DWDZ(c,t)*n2*n3*a23+
					C_DWDX(c,t)*n3*n3*a31+C_DWDY(c,t)*n3*n3*a32+C_DWDZ(c,t)*n3*n3*a33)
							
						)*f;


phi_33_w2=phi_33_w2+
			c2bs*C_K(c,t)*(
			
			C_DUDX(c,t)*n1*n1+C_DUDY(c,t)*n1*n2+C_DUDZ(c,t)*n1*n3+
			C_DVDX(c,t)*n2*n1+C_DVDY(c,t)*n2*n2+C_DVDZ(c,t)*n2*n3+
			C_DWDX(c,t)*n3*n1+C_DWDY(c,t)*n3*n2+C_DWDZ(c,t)*n3*n3
		
						)

			*(n3*n3-1./3.*(n1*n1+n2*n2+n3*n3))*f;


// prüfen ob die Dichte berücksichtigt werden muss!!

	dS[eqn]= -c1s*C_D(c,t)/C_K(c,t)*2.*n3*n3*f;
	source= phi_33_w1+phi_33_w2;

return source;
}


DEFINE_SOURCE(uv_source, c, t, dS, eqn)
{
real source;
real u11, u12, u13, u21, u22, u23, u31, u32, u33;
real a11, a12, a13, a21, a22, a23, a31, a32, a33;
real n1, n2, n3, f;
real phi_12_w1, phi_12_w2;

f=pow(C_K(c,t),1.5)/C_D(c,t)/cl/C_UDSI(c,t,0);


n1=C_UDMI(c,t,0);
n2=C_UDMI(c,t,1);
n3=C_UDMI(c,t,2);


u11=C_RUU(c,t);
u12=C_RUV(c,t);
u13=C_RUW(c,t);
u21=C_RUV(c,t);
u22=C_RVV(c,t);
u23=C_RVW(c,t);
u31=C_RUW(c,t);
u32=C_RVW(c,t);
u33=C_RWW(c,t);

a11=u11/C_K(c,t)-2./3.;
a12=u12/C_K(c,t);
a13=u13/C_K(c,t);
a21=u21/C_K(c,t);
a22=u22/C_K(c,t)-2./3.;
a23=u23/C_K(c,t);
a31=u31/C_K(c,t);
a32=u32/C_K(c,t);
a33=u33/C_K(c,t)-2./3.;



phi_12_w1=c1s*C_D(c,t)/C_K(c,t)*(
								 -3./2.*(u11*n1*n2+u21*n2*n2+u31*n3*n2)
								 -3./2.*(u12*n1*n1+u22*n1*n2+u32*n1*n3)
								 )*f;
phi_12_w2=-c2s*(
				C_DUDX(c,t)*u11+C_DVDX(c,t)*u21+C_DWDX(c,t)*u31+
				C_DUDY(c,t)*u12+C_DVDY(c,t)*u22+C_DWDY(c,t)*u32+
				C_DUDZ(c,t)*u13+C_DVDZ(c,t)*u23+C_DWDZ(c,t)*u33
				)*(-3.*n1*n2)*f;

phi_12_w2=phi_12_w2-
			c2as*C_K(c,t)*(

			-3./2.*(C_DUDX(c,t)*n1*n2*a11+C_DUDY(c,t)*n1*n2*a12+C_DUDZ(c,t)*n1*n2*a13+
					C_DUDX(c,t)*n2*n2*a21+C_DUDY(c,t)*n2*n2*a22+C_DUDZ(c,t)*n2*n2*a23+
					C_DUDX(c,t)*n3*n2*a31+C_DUDY(c,t)*n3*n2*a32+C_DUDZ(c,t)*n3*n2*a33)

			-3./2.*(C_DVDX(c,t)*n1*n1*a11+C_DVDY(c,t)*n1*n1*a12+C_DVDZ(c,t)*n1*n1*a13+
					C_DVDX(c,t)*n2*n1*a21+C_DVDY(c,t)*n2*n1*a22+C_DVDZ(c,t)*n2*n1*a23+
					C_DVDX(c,t)*n3*n1*a31+C_DVDY(c,t)*n3*n1*a32+C_DVDZ(c,t)*n3*n1*a33)
							
						)*f;


phi_12_w2=phi_12_w2+
			c2bs*C_K(c,t)*(
			
			C_DUDX(c,t)*n1*n1+C_DUDY(c,t)*n1*n2+C_DUDZ(c,t)*n1*n3+
			C_DVDX(c,t)*n2*n1+C_DVDY(c,t)*n2*n2+C_DVDZ(c,t)*n2*n3+
			C_DWDX(c,t)*n3*n1+C_DWDY(c,t)*n3*n2+C_DWDZ(c,t)*n3*n3
		
						)

			*(n1*n2)*f;


// prüfen ob die Dichte berücksichtigt werden muss!!

	dS[eqn]= -c1s*C_D(c,t)/C_K(c,t)*(3./2.*n1*n1+3./2.*n2*n2)*f;
	source= phi_12_w1+phi_12_w2;

return source;
}

DEFINE_SOURCE(uw_source, c, t, dS, eqn)
{
real source;
real u11, u12, u13, u21, u22, u23, u31, u32, u33;
real a11, a12, a13, a21, a22, a23, a31, a32, a33;
real n1, n2, n3, f;
real phi_13_w1, phi_13_w2;

f=pow(C_K(c,t),1.5)/C_D(c,t)/cl/C_UDSI(c,t,0);


n1=C_UDMI(c,t,0);
n2=C_UDMI(c,t,1);
n3=C_UDMI(c,t,2);


u11=C_RUU(c,t);
u12=C_RUV(c,t);
u13=C_RUW(c,t);
u21=C_RUV(c,t);
u22=C_RVV(c,t);
u23=C_RVW(c,t);
u31=C_RUW(c,t);
u32=C_RVW(c,t);
u33=C_RWW(c,t);

a11=u11/C_K(c,t)-2./3.;
a12=u12/C_K(c,t);
a13=u13/C_K(c,t);
a21=u21/C_K(c,t);
a22=u22/C_K(c,t)-2./3.;
a23=u23/C_K(c,t);
a31=u31/C_K(c,t);
a32=u32/C_K(c,t);
a33=u33/C_K(c,t)-2./3.;



phi_13_w1=c1s*C_D(c,t)/C_K(c,t)*(
								 -3./2.*(u11*n1*n3+u21*n2*n3+u31*n3*n3)
								 -3./2.*(u13*n1*n1+u23*n1*n2+u33*n1*n3)
								 )*f;
phi_13_w2=-c2s*(
				C_DUDX(c,t)*u11+C_DVDX(c,t)*u21+C_DWDX(c,t)*u31+
				C_DUDY(c,t)*u12+C_DVDY(c,t)*u22+C_DWDY(c,t)*u32+
				C_DUDZ(c,t)*u13+C_DVDZ(c,t)*u23+C_DWDZ(c,t)*u33
				)*(-3.*n1*n3)*f;

phi_13_w2=phi_13_w2-
			c2as*C_K(c,t)*(

			-3./2.*(C_DUDX(c,t)*n1*n3*a11+C_DUDY(c,t)*n1*n3*a12+C_DUDZ(c,t)*n1*n3*a13+
					C_DUDX(c,t)*n2*n3*a21+C_DUDY(c,t)*n2*n3*a22+C_DUDZ(c,t)*n2*n3*a23+
					C_DUDX(c,t)*n3*n3*a31+C_DUDY(c,t)*n3*n3*a32+C_DUDZ(c,t)*n3*n3*a33)

			-3./2.*(C_DWDX(c,t)*n1*n1*a11+C_DWDY(c,t)*n1*n1*a12+C_DWDZ(c,t)*n1*n1*a13+
					C_DWDX(c,t)*n2*n1*a21+C_DWDY(c,t)*n2*n1*a22+C_DWDZ(c,t)*n2*n1*a23+
					C_DWDX(c,t)*n3*n1*a31+C_DWDY(c,t)*n3*n1*a32+C_DWDZ(c,t)*n3*n1*a33)
							
						)*f;


phi_13_w2=phi_13_w2+
			c2bs*C_K(c,t)*(
			
			C_DUDX(c,t)*n1*n1+C_DUDY(c,t)*n1*n2+C_DUDZ(c,t)*n1*n3+
			C_DVDX(c,t)*n2*n1+C_DVDY(c,t)*n2*n2+C_DVDZ(c,t)*n2*n3+
			C_DWDX(c,t)*n3*n1+C_DWDY(c,t)*n3*n2+C_DWDZ(c,t)*n3*n3
		
						)

			*(n1*n3)*f;


// prüfen ob die Dichte berücksichtigt werden muss!!

	dS[eqn]= -c1s*C_D(c,t)/C_K(c,t)*(3./2.*n1*n1+3./2.*n3*n3)*f;
	source= phi_13_w1+phi_13_w2;

return source;
}


DEFINE_SOURCE(vw_source, c, t, dS, eqn)
{
real source;
real u11, u12, u13, u21, u22, u23, u31, u32, u33;
real a11, a12, a13, a21, a22, a23, a31, a32, a33;
real n1, n2, n3, f;
real phi_23_w1, phi_23_w2;

f=pow(C_K(c,t),1.5)/C_D(c,t)/cl/C_UDSI(c,t,0);


n1=C_UDMI(c,t,0);
n2=C_UDMI(c,t,1);
n3=C_UDMI(c,t,2);


u11=C_RUU(c,t);
u12=C_RUV(c,t);
u13=C_RUW(c,t);
u21=C_RUV(c,t);
u22=C_RVV(c,t);
u23=C_RVW(c,t);
u31=C_RUW(c,t);
u32=C_RVW(c,t);
u33=C_RWW(c,t);

a11=u11/C_K(c,t)-2./3.;
a12=u12/C_K(c,t);
a13=u13/C_K(c,t);
a21=u21/C_K(c,t);
a22=u22/C_K(c,t)-2./3.;
a23=u23/C_K(c,t);
a31=u31/C_K(c,t);
a32=u32/C_K(c,t);
a33=u33/C_K(c,t)-2./3.;



phi_23_w1=c1s*C_D(c,t)/C_K(c,t)*(
								 -3./2.*(u12*n1*n3+u22*n2*n3+u32*n3*n3)
								 -3./2.*(u13*n2*n1+u23*n2*n2+u33*n2*n3)
								 )*f;
phi_23_w2=-c2s*(
				C_DUDX(c,t)*u11+C_DVDX(c,t)*u21+C_DWDX(c,t)*u31+
				C_DUDY(c,t)*u12+C_DVDY(c,t)*u22+C_DWDY(c,t)*u32+
				C_DUDZ(c,t)*u13+C_DVDZ(c,t)*u23+C_DWDZ(c,t)*u33
				)*(-3.*n2*n3)*f;

phi_23_w2=phi_23_w2-
			c2as*C_K(c,t)*(

			-3./2.*(C_DVDX(c,t)*n1*n3*a11+C_DVDY(c,t)*n1*n3*a12+C_DVDZ(c,t)*n1*n3*a13+
					C_DVDX(c,t)*n2*n3*a21+C_DVDY(c,t)*n2*n3*a22+C_DVDZ(c,t)*n2*n3*a23+
					C_DVDX(c,t)*n3*n3*a31+C_DVDY(c,t)*n3*n3*a32+C_DVDZ(c,t)*n3*n3*a33)

			-3./2.*(C_DWDX(c,t)*n1*n2*a11+C_DWDY(c,t)*n1*n2*a12+C_DWDZ(c,t)*n1*n2*a13+
					C_DWDX(c,t)*n2*n2*a21+C_DWDY(c,t)*n2*n2*a22+C_DWDZ(c,t)*n2*n2*a23+
					C_DWDX(c,t)*n3*n2*a31+C_DWDY(c,t)*n3*n2*a32+C_DWDZ(c,t)*n3*n2*a33)
							
						)*f;


phi_23_w2=phi_23_w2+
			c2bs*C_K(c,t)*(
			
			C_DUDX(c,t)*n1*n1+C_DUDY(c,t)*n1*n2+C_DUDZ(c,t)*n1*n3+
			C_DVDX(c,t)*n2*n1+C_DVDY(c,t)*n2*n2+C_DVDZ(c,t)*n2*n3+
			C_DWDX(c,t)*n3*n1+C_DWDY(c,t)*n3*n2+C_DWDZ(c,t)*n3*n3
		
						)

			*(n2*n3)*f;


// prüfen ob die Dichte berücksichtigt werden muss!!

	dS[eqn]= -c1s*C_D(c,t)/C_K(c,t)*(3./2.*n2*n2+3./2.*n3*n3)*f;
	source= phi_23_w1+phi_23_w2;

return source;
}


DEFINE_ON_DEMAND(wall_dist_demand)
{
Domain *d; /* declare domain pointer since it is not passed as an argument to the DEFINE macro */

Thread *t;
cell_t c;

d = Get_Domain(1); /* Get the domain using Fluent utility */

thread_loop_c(t,d)
{
begin_c_loop(c,t)
{
				C_UDSI(c,t,0) = C_WALL_DIST(c,t);
}
end_c_loop(c,t)
}
}


DEFINE_ON_DEMAND(n_demand)
{
	real n1, n2, n3, n[ND_ND];


Domain *d; /* declare domain pointer since it is not passed as an argument to the DEFINE macro */

Thread *t;
cell_t c;

d = Get_Domain(1); /* Get the domain using Fluent utility */

thread_loop_c(t,d)
{
begin_c_loop(c,t)
{
				NV_V(n,=,C_UDSI_G(c,t,0));
				C_UDMI(c,t,0)=n[0]/NV_MAG(n);
				C_UDMI(c,t,1)=n[1]/NV_MAG(n);
				C_UDMI(c,t,2)=n[2]/NV_MAG(n);
}
end_c_loop(c,t)
}
}