Back to report
#include#define m 1.5 #define d (0.8724*(m+1)) #define R_e 43.48 #define P_r 1.8 #define R_m 404 #define B_0 0.1827 #define a (3.14159265/d) /*"2pi/2d" : Alpha */ void difeq(int k, int k1, int k2, int jsf, int is1, int isf, int indexv[], int ne, float **s, float **y) { /*initializes all **s to 0*/ int i1,i2; for (i1 = 1; i1 < 10; i1++) for (i2 = 1; i2 < 19; i2++) s[i1][i2]=0.0; for (i1 = 1; i1 < 10; i1++) s[i1][jsf]=0.0; if (k==k1) { s[5][10]=1.0; s[5][jsf]=y[1][1]; s[6][11]=1.0; s[6][jsf]=y[2][1]; s[7][12]=1.0; s[7][jsf]=y[3][1]; s[8][13]=1.0; s[8][jsf]=y[4][1]; s[9][14]=1.0; s[9][jsf]=y[5][1] + 1.0; } else if (k<=k2) { /*"h=d/M"*/ float h=d/(k2-1); /*"z[k]=(k-1)*d/M"*/ float z_k=(k-1)*h; /*"z[k-1]=(k-2)*d/M"*/ float z_k1=(k-2)*h; /* val = temperature at z_k = pressure scale height at z_k */ float val=(m+1.0-z_k)/(m+1.0); /* val1 = temperature at z_k1 = pressure scale height at z_k1 */ float val1=(m+1.0-z_k1)/(m+1.0); /* avg "-1/val1 & -1/val" */ float dS_0=(-0.5)*((1.0/val1)+(1.0/val)); /* (-m/(m+1)) * avg "-1/val1 & -1/val" */ float ovh_0=(-m/(m+1.0))*dS_0; /* (1/(m+1)) * avg "-1/val1 & -1/val" */ float dT_0=(1.0/(m+1.0))*dS_0; /*avg"val1^m & val^m"*/ float rho_0=(0.5)*(pow(val1,m)+pow(val,m)); /*avg"val1^-m & val^-m"*/ float ovrho_0=(0.5)*(pow(val1,-m)+pow(val,-m)); /*avg y[1] valueX2*/ float Y1=(y[1][k]+y[1][k-1]); /*avg y[2] valueX2*/ float Y2=(y[2][k]+y[2][k-1]) ; /*avg y[3] valueX2*/ float Y3=(y[3][k]+y[3][k-1]); /*avg y[4] valueX2*/ float Y4=(y[4][k]+y[4][k-1]); /*avg y[5] valueX2*/ float Y5=(y[5][k]+y[5][k-1]); /*avg y[6] valueX2*/ float Y6=(y[6][k]+y[6][k-1]); /*avg y[7] valueX2*/ float Y7=(y[7][k]+y[7][k-1]); /*avg y[8] valueX2*/ float Y8=(y[8][k]+y[8][k-1]); /*avg y[9] value*/ float Y9=0.5*(y[9][k]+y[9][k-1]); /*Gamma*/ float Y=5.0/3.0; s[1][1]=-1.0; s[1][5]=-0.5*h; s[1][10]=1.0; s[1][14]=-0.5*h; s[1][jsf]=y[1][k]-y[1][k-1]-h*0.5*Y5; s[2][2]=-1.0; s[2][6]=-0.5*h; s[2][11]=1.0; s[2][15]=-0.5*h; s[2][jsf]=y[2][k]-y[2][k-1]-h*0.5*Y6; s[3][3]=-1.0; s[3][7]=-0.5*h; s[3][12]=1.0; s[3][16]=-0.5*h; s[3][jsf]=y[3][k]-y[3][k-1]-h*0.5*Y7; s[4][4]=-1.0; s[4][8]=-0.5*h; s[4][13]=1.0; s[4][17]=-0.5*h; s[4][jsf]=y[4][k]-y[4][k-1]-h*0.5*Y8; s[5][1]=-0.5*h*(a*a+R_e*rho_0*Y9); s[5][2]=0.5*h*((a*a*ovrho_0)*(R_e*R_m*B_0*B_0+(4.0/3.0)*ovh_0*ovh_0)); s[5][3]=0.5*h*(R_e*rho_0*a); s[5][4]=-0.5*h*(R_e*R_m*a*Y9*B_0); s[5][5]=-1.0+0.5*h*ovh_0; s[5][9]=-0.25*h*(R_e*rho_0*Y1); s[5][10]=s[5][1]; s[5][11]=s[5][2]; s[5][12]=s[5][3]; s[5][13]=s[5][4]; s[5][14]=1.0+0.5*h*ovh_0; s[5][18]=s[5][9]; s[5][jsf]=(y[5][k]-y[5][k-1])+(Y1*s[5][1])+(Y2*s[5][2])+(Y3*s[5][3])+(Y4*s[5][4])+(0.5*h*ovh_0)*(Y5); s[6][1]=-0.5*h*rho_0; s[6][2]=-0.5*h*a*a; s[6][6]=-1.0+0.5*h*ovh_0; s[6][10]=s[6][1]; s[6][11]=s[6][2]; s[6][15]=1.0+0.5*h*ovh_0; s[6][jsf]=(y[6][k]-y[6][k-1])+(Y1*s[6][1])+(Y2*s[6][2])+(0.5*h*ovh_0)*(Y6); s[7][2]=-0.5*h*(P_r*R_e*a)*(-dS_0+(0.5*Y*B_0*B_0)*ovrho_0*ovh_0*ovh_0); s[7][3]=-0.5*h*(a*a+P_r*R_e*Y9*rho_0); s[7][7]=-1.0+(0.5*h)*dT_0; s[7][9]=-0.25*h*(P_r*R_e*rho_0*Y3); s[7][11]=s[7][2]; s[7][12]=s[7][3]; s[7][16]=1.0+(0.5*h)*dT_0; s[7][18]=s[7][9]; s[7][jsf]=(y[7][k]-y[7][k-1])+(Y2*s[7][2])+(Y3*s[7][3])+(0.5*h)*dT_0*Y7; s[8][2]=0.5*h*(a*R_m*B_0*ovrho_0); s[8][4]=-0.5*h*(a*a+R_m*Y9); s[8][8]=-1.0; s[8][9]=-0.25*h*(R_m*Y4); s[8][11]=s[8][2]; s[8][13]=s[8][4]; s[8][17]=1.0; s[8][18]=s[8][9]; s[8][jsf]=(y[8][k]-y[8][k-1])+(Y2*s[8][2])+(Y4*s[8][4]); s[9][9]=-1.0; s[9][18]=1.0; s[9][jsf]=y[9][k]-y[9][k-1]; } else { s[1][10]=1.0; s[1][jsf]=y[1][k2]; s[2][11]=1.0; s[2][jsf]=y[2][k2]; s[3][12]=1.0; s[3][jsf]=y[3][k2]; s[4][13]=1.0; s[4][jsf]=y[4][k2]; } }
This code is based on an algorithm detailed in Numerical Recipes in c in chapter 16
Numerical Recipes in c by Press, Flannery, Teukolsky, Vetterling
Cambridge University Press, 1988