unsteady heat flux

co2

following is the udf that I am using to define heat flux at a surface in my problem. Neglect the different equations, but can you please tell me if my udf is syntactically correct? i am getting very weird results and impossibly high velocities .. i was just wondering if someone can spot an obvious stupid mistake in the following code before I dig into other possibilities .. This is just my second time writing a UDF and hence I am not experienced at all. Hence requesting your help .. thanks ..

#include "udf.h"

DEFINE_PROFILE(unsteady_heatflux, thread, position) {

face_t f;

begin_f_loop(f, thread)

{

int n; // it is the day of the year.

// all the important angles in terms of degrees.

float phi, delta, beta, gamma, omega;

// constants for solar time calculation.

float E,B,solar_t;

float heat_flux[24];

char character;

int i;

float value;

float G_Constant, G_ans;

float cos_theta;

float longitude_std = 120.0; // this is the longitude in degrees for the standard Pacific time where Las Vegas is in.

float longitude_local = 115.0+(9.0*1.0/60.0)+(19.0*1.0/3600.0); // this is longitude in degrees for las vegas where the station is.

float pi = 3.141592654;

int hour;

real t = RP_Get_Real("flow-time");

/*************************/

beta=0.0;

beta=beta*(pi/180.0);

gamma=0.0;

gamma=gamma*(pi/180.0);

/************************/

phi = 36.0+(4.0*1.0/60.0)+(44.0*1.0/3600.0) ;

phi=phi*(pi/180.0);

heat_flux[0]=0.0000;

heat_flux[1]=0.0000;

heat_flux[2]=0.0000;

heat_flux[3]=0.0000;

heat_flux[4]=0.0000;

heat_flux[5]=0.0000;

heat_flux[6]=0.0000;

heat_flux[7]=1.1000;

heat_flux[8]=90.4100;

heat_flux[9]=239.3100;

heat_flux[10]=397.5100;

heat_flux[11]=467.2100;

heat_flux[12]=394.8100;

heat_flux[13]=417.9100;

heat_flux[14]=256.7000;

heat_flux[15]=190.3000;

heat_flux[16]=80.8000;

heat_flux[17]=3.6000;

heat_flux[18]=0.0000;

heat_flux[19]=0.0000;

heat_flux[20]=0.0000;

heat_flux[21]=0.0000;

heat_flux[22]=0.0000;

heat_flux[23]=0.0000;

n=(t/86400)+1; // this is a simple equation to find the day of the year from t. 86400 sec in one day.

B=(360.0*(n-81)/364.0)*(pi/180.0); // B is basically some angle but since c language requires angles to be in radian, we are

// converting it here into radians.

E = 9.87*sin(2.0*B) - 7.53*cos(B) - 1.5*sin(B);

solar_t = t + 4.0*(longitude_std-longitude_local)*60.0+E*60.0; // notice that we multiply by 60 to get everything in seconds.

if((solar_t-((n-1)*86400))<=43200)

omega = -((43200.0-solar_t)/3600.0)*15.0;

else

omega = ((solar_t)/3600.0)*15.0;

delta = 23.45*sin(360.0*((284+n)/365.0)*(pi/180.0));

// it is now necessary to find out which hour of the day we are in at a particular point during the solution.

// 1 hr = 3600 sec.

hour = ((t-((n-1)*86400))/3600.0) + 1; // will give us the integer part of the division. thus in the first 3600 seconds of the day we are in the

// first hour.

omega=omega*(pi/180.0);

delta=delta*(pi/180.0);

G_Constant = heat_flux[hour-1]/(sin(delta)*sin(phi) + cos(delta)*cos(phi)*cos(omega));

cos_theta=sin(delta)*sin(phi)*cos(beta) - sin(delta)*cos(phi)*sin(beta)*cos(gamma) + cos(delta)*cos(phi)*cos(beta)*cos(omega) + cos(delta)*sin(phi)*sin(beta)*cos(gamma)*cos(omega ) + cos(delta)*sin(beta)*sin(gamma)*sin(omega);

G_ans= G_Constant*cos_theta;

if(G_ans<0.0)

G_ans=0.0;

F_PROFILE(f, thread, position) = G_ans;

} end_f_loop(f, thread) }