laserlengths.cpp

this progam calculates the optimum laser powers at different rod lengths

/************************************
this progam calculates the optimum laser powers at different rod lengths. the user input a rod length, the program
 generates corresponding laser power. then input another rod length, run the program again to get its laser power. 
 finally, the user can get teh relationship, optimum laser power over rod length. the user input Lby setting L(at 
 the sixth line of the program)to the desired value. also the N(at the fourth line) is reset to: N=500*L for different L.
 */
#include<math.h>
#include<iomanip.h>
#include<fstream.h>

const int N=300, M=10, Iter_Num=4;
const long double L=0.6, R=0.04, tip_r=0.015, hight=0.02, a=0.002, Ea=182004.0, Rgas=8.314, K0=237000.0, delta_t= 0.002, w=0.01;
void SolveTrid(long double [N+1], long double [N+1], long double [N+1], long double [N+1], long double [N+1]);
void TempDistr(long double, long double, long double, long double[N+1]);
void GetTgrowth(long double [M+1]);
long double growth (long double, long double, long double, long double);
void One_iterat (long double, long double &);

void main()
{
	long double P0[Iter_Num];
	P0[0]=1.48;  //initial guess for P

	for (int i=0; i<Iter_Num; i++);
	{
		One_iterat(P0[i],P0[i+1]);
    	cout<<"After"<<i+1<<"iteration(s),P0 is:"<<setiosflags(ios::showpoint)<<setprecision(20)<<P0[i+1]<<endl;
      
	}
}
/**************************************************************
this following subroutine is going through one cycle from step 1 to step 7 in section 5.6 of chapter five,
calculating Pk from Pk-1
*****************************************************************/
void One_iterat(long double P0, long double &new_P)
{
	long double T[M+1]={0}, Pdelta_T[M+1]={0}, X[M+1]={0},XXtX[M+1], T_Tgrow[M+1], Tgrow[M+1];
	long double delta_P0, temp1=0.0, XXt=0.0, S=0.0, Uk=2;

	delta_P0=P0/1000;

	GetTgrowth(Tgrow);

	TempDistr(L,P0,w,T);

	TempDistr(L, P0+delta_P0, w, Pdelta_T);
	for(int n=0; n<M+1; n++)
		X[n]=(Pdelta_T[n]-T[n])/delta_P0;

	for(int i=0;i<M+1; i++)
		XXt=XXt+X[i]*X[i];
	for(i=0; i<M+1;i++)
		XXtX[i]=X[i]/XXt;

	for(i=0; i<M+1;i++)
	{
		T_Tgrow[i]=T[i]-Tgrow[i];
		S=S+(T_Tgrow[i]*T_Tgrow[i]);
	}
	cout<<endl<<"the least square norm at this iteration is:"<<S<<endl;
	temp1=0.0;
	for(i=0;i<M+1;i++)
		temp1=temp1+XXtX[i]*T_Tgrow[i];
	new_P=P0+temp1;
}

/**********************************************************
the following subroutine calculates the temperature distribtuion on the whole surface of the rod
******************************************************/
void TempDistr(long double L,long double P0, long double w, long double Temp[M+1])
{
	long double L1, c, h, Lumda;
	long double Ks=0.0017, absor =1.0, f, kgas=0.001, Kd=1.65, Nu=0.36, pai=3.1415926;
	long double x[N+1], r[N+1], O[N+1], x_half[N+1], sq_r_half[N+1], h_conv[N+1], Qin[N+1], b[N+1], slope[N+1];
	long double U[N+1]={0}, l[N+1]={0}, d[N+1]={0};
    long double y=R*R-tip_r*tip_r;
	L1=L-hight;
	c=sqrt(y/hight);
	h=L/N;
	Lumda=2*absor/(pai*w*w);

	cout<<endl;
	for(int i=0; i<N+1; i++)
	{
		x[i]=i*h;
		if(x[i]<=L1)
			r[i]=R;
		else
		{
			r[i]=c*sqrt(L+tip_r*tip_r/(c*c)-x[i]);
			slope[i]= c*c/(2*r[i]);
		}
		h_conv[i]=Nu*kgas/(2*r[i]);
	}
	
	for(i=1; i<N+1;i++)
	{
		x_half[i]=i*h-h/2;
		if(x_half[i]<=L1)
			sq_r_half[i]=R*R;
		else
			sq_r_half[i]=c*sqrt(L+tip_r/(c*c)-x_half[i])*c*sqrt(L+tip_r*tip_r/(c*c)-x_half[i]);
	}

	d[0]=3*Ks/(8*Kd*r[0]*r[0]*r[0])+1/h;
	U[0]=-1/h;
	l[N]=-1/h;
	d[N]=h_conv[N]/(pai*Kd*tip_r*tip_r)+1/h;

	for(i=1;i<N; i++)
	{
		l[i]=-(Kd*sq_r_half[i]/h)/h;
		d[i]=-r[i]*h_conv[i]+Kd*(sq_r_half[i]+sq_r_half[i+1])/(h*h);
		U[i]=-(Kd*sq_r_half[i+1]/h)/h;
	}

	for(i=1;i<N;i++)
	{
		if(x[i]<=L1)
			Qin[i]=0.0;
		else
			Qin[i]=2*P0*absor*(x[i]/sqrt(slope[i]+1))*exp(-2*r[i]*r[i]/(w*w))/(pai*w*w);
	}
	
	b[0]=0.0;
	b[N]=Lumda*P0*exp(-2*tip_r*tip_r/(w*w))/(pai*Kd*tip_r*tip_r);

	for(i=1;i<N;i++)
		b[i]=r[i]*Qin[i];
	
	SolveTrid(O,b,l,U,d);

	for(i=0;i<=N+1;i++)
		O[i]=O[i]+300;
	
	for(i=0;i<M+1;i++)
		Temp[i]=O[N+i-M];
}

/***********************************************************************
the following subroutine solves the linear system by Thomas Algorithm
************************************************/

void SolveTrid(long double O[N+1], long double b[N+1], long double l[N+1], long double U[N+1], long double d[N+1])
{
	long double b1[N+1],d1[N+1];
	d1[0]=d[0];
	b1[0]=b[0];
	for(int i=1; i<N+1; i++)
	{
		d1[i]= d[i]-l[i]*U[i-1]/d1[i-1];
		b1[i]= b[i]-l[i]*b1[i-1]/d1[i-1];
	}
	O[N]=b1[N]/d1[N];

	for(i=N-1;i>=0; i--)
		O[i]=(b1[i]-U[i]*O[i+1])/d1[i];
}

/************************************************************************
the following subroutine calculates the required temperature distribution Td onthe parabolic portion of the rod.
**********************************************************************/

void GetTgrowth(long double Td[M+1])
{
	long double L1,e,f,c,h;
	long double x[M+1],r[M+1],slope[M+1],x_new1[M+1], r_new1[M+1], x_new2[M+1], r_new2[M+1],delta_x[M+1];

	long double y=R*R-tip_r*tip_r;

	L1=L-hight;
	c=sqrt(y/hight);
	h=hight/M;

	cout<<endl;
	for(int i=0; i<M+1; i++)
	{
		x[i]=L1+i*h;
		r[i]=c*sqrt(L+tip_r*tip_r/(c*c)-x[i]);
		slope[i]=c*c/(2*r[i]);
	}
	delta_x[0]=0;
	delta_x[M]=a;

	for(int j=1; j<M; j++)
	{
		r_new1[j]=R;
		x_new1[j]=slope[j]*(R-r[j])+x[j];
		if(x_new1[j]<=(L1+a))
			delta_x[j]=growth(r[j], x[j], r_new1[j], x_new1[j]);
		else break;
	}

	for(int k=j; k<M; k++)
	{
		e=(x[k]-slope[k]*r[k]-L-a-tip_r*tip_r/(c*c))/(c*c);
		f=sqrt(slope[k]-4*e);
		r_new2[k]=c*c*(f-slope[k])/2;
		x_new2[k]=slope[k]*(r_new2[k]-r[k])+x[k];
		delta_x[k]=growth(r[k],x[k],r_new2[k], x_new2[k]);
	}

	for(int p=0;p<M+1;p++)
		Td[p]=-Ea/(Rgas*log(delta_x[p]/(K0*delta_t)));
	Td[0]=Td[1]-(Td[2]-Td[1])*(Td[2]-Td[1])/(Td[3]-Td[2]);
}
/***************************************************
the following subroutine calculates the growth of each grid point from current layer to the next layer
*****************************************************/
long double growth(long double a1, long double b1, long double a2, long double b2)
{
	long double d;
	d=(a2-a1)*(a2-a1)+(b2-b1)*(b2-b1);
	return sqrt(d);
}