usfkad6_05.cpp

Solutions are obtained for Poissson, diffusion, or wave PDEs homogeneous or nonhomogeneous equations

	{
		term = term+1;
	}
	if((mtrx[t][6]=="N"))  //(t != 4)&& omitted
	{
		term = term+1;
	}
   } 

//   In case there's no steady state term.
   if (term==0)
   {
	   outFilE << 0 << "\\\\ ";
   }
   else
   {
   for(i=1; i<term; i++)
   {
	   itoa (i,sterm, 10);
	   outFilE << "\\Psi_" << i <<  "+";
   }

itoa (term,sterm, 10);
outFilE << "\\Psi_" << term << "\\\\ ";
   }

if(eqtype[0]==3)
{
	outFilE << "\\Psi_{transient}=\\Psi_{transient \\, \\#1}+\\Psi_{transient \\, \\#2}" << "\\\\";
}
else
{
	outFilE << "\\\\" ;
}

// Finished first line of output *******************************************
   
/* set the number of terms counter to 0 */
term=0;

/* Look for low end, high end nonhomogeneities*/
for(lowhigh=1; lowhigh<3; lowhigh++)
       {
	   if(lowhigh==1)  // low end
	     {
           lohiint = 2;
		   lohisym = "NH";
		   p=8;
	   }
		if(lowhigh==2)  // high end
	     {
           lohiint = 6;
		   lohisym = "HN";
		   p = 9;
		}

for(i=1;i<dim+1;i++) // i loops through dim
 {
  /* identify the coordinate row r(i)*/
  r=eqtype[2+i];
  if (mtrx[r][lohiint]=="N")  //Is row(i) nonhomog?
  {
	  term=term+1; 
      itoa(term, sterm,10);
	  outFilE << "\n \\\\ \\\\ \\\\ \\\\ \\Psi_" << term << "=";
	  
// Re-initialize the output strings
Symb1="";
Symb2="=";
Eigans1="";
Nonans1="\\; ";
Eigans2="\n \\\\ where";
// Nonans2="\n\\\\ ";
Nonans2="\\\\ ";
Eigans53="";
Nonans3="\\; ";
Eigans4="\n\\\\ ";
Nonans4="\\\\ ";

  // ith coord, in row r, is nonhomog; Call the eigenfunctions
  for (k=1;k<dim+1;k++)
    {
     if (k!=i) // k loops through the homog rows
       {
        re=eqtype[k+2]; 
        // re(k) is the row number for the current eigenfunction
//   Change to Caucy-Euler equation if 2D and steadystate
/*		if (re==5 && eqtype[1]==2 && (eqtype[0]==0 || eqtype[0] == 1 || eqtype[0] == 3))
		{
			re=4;
		}                   */
		
		f4 = mtrx[re][5];
		if (f4=="I" || f4=="O")
		{
			f4="S";
		}
		
		f=mtrx[re][3] + mtrx[re][1] + f4 + "HH";

        A.ans1 = f;
        A.ans2 = mtrx[re][7];
        A.ans3 = mtrx[re][9];
        A.symb1 = "" ;
        A.symb2 = "" ; 
        /* assemble the answer for this condition */

		//   Debug
/*		  cout<<f<<"\n";
  cout<<r<<"\n";
  cout<<lohisym<<"\n";
*/

		A=sul(A,r);
		
		Symb1 += A.symb1;
		Symb2 += A.symb2;
		Eigans1 += A.ans1;
		Eigans2 += "\\\\ \n" + A.ans2;
		Eigans53 += "\\;" + A.ans5 + "\\;" + A.ans3;
		Eigans4 += A.ans4 + "\\\\";
	// That builds the answer string for this factor.
       }
     }
  // Insert the nonhomogeneous factor for variable i

//   Change to Caucy-Euler equation if 2D and steadystate
/*		if (r==5 && eqtype[1]==2 && (eqtype[0]==0 || eqtype[0] == 1 || eqtype[0] == 3))
		{
			r=4;
		}              */
		
		A=KandF(A,mtrx[r][p], r); // (A, nonhomo variable, its symbol); KandF packs KAPPA, AAA, & nonhomfcn
/*		if (mtrx[r][5]=="O")  // Modify for outgoing.
		{
			A=KandF(A,"O", r);
		}
		if (mtrx[r][5]=="I")  // Modify for incoming.
		{
			A=KandF(A,"I", r);
		}
*/
  
 nf=mtrx[r][3] + mtrx[r][1]+ mtrx[r][5] + lohisym;
  A.ans1 = nf;
  A.ans2= mtrx[r][7];
  A.ans3 = mtrx[r][9];

  // Debug
/*  cout<<nf<<"\n";
  cout<<r<<"\n";
  cout<<lohisym<<"\n";  */
  A = sul(A,r);
  
//  Add the nonhomogeneous factor strings

Nonans1+=A.ans1;
Nonans2 += A.ans2;
Nonans3 += A.ans3;
Nonans4+=A.ans4;
 
//   Write to outfile
outFilE << Symb1 << Eigans1 << Nonans1 << "\\;" << A.AAA << Eigans2 << Nonans2 << "\\\\";
outFilE << A.AAA << Symb2 << Eigans53 <<  Nonans3 << "\\;" << A.nonhomf << Eigans4 << Nonans4;
  }  // End "i nonhomo"
}   // End i = 1 to dim
        }  //  End lowhigh loop

//    Clear "r"
/*  r=100;  */

//  Append Green's function term
if(eqtype[7]==1 || eqtype[0]==2 || eqtype[0]==4)
{
	A=Green(A,term);
// Green writes to outfile
}

// Append transient terms
if(eqtype[0]==1)  // Diffusion
{
	t=0;
	A=Trans(A,t);
}
if(eqtype[0]==3)  //  Wave
{
	t=1;
	A=Trans(A,t);
}
// Trans writes to outfile

outFilE << "$\n\\end{document}";
outFilE.close();
cout << "DONE\n";
}


	pdeans Green(pdeans A, int term)
//pdeans Green(pdeans A, std::string sterm)
// Assembles and writes Green's function term to the outfile
{
std::string KAPPA, nonhomf, AAA, f, f4;
std::string Symb1, Symb2, Eigans1, Nonans1, Eigans2, Nonans2, Eigans53, Nonans3, Eigans4, Nonans4;
term = term+1;
outFilE << "\n  \\\\ \\\\ \\\\ \\Psi_{" << term << "} = ";
int k, re;

// Re-initialize the output strings
Symb1="";
Symb2="=";
Eigans1="";
Eigans2="\n\\\\where";
Eigans53="";
Eigans4="\n\\\\";

  for (k=1;k<eqtype[1]+1;k++)
    {
        re=eqtype[k+2]; 
        // re(k) is the row number for the current eigenfunction
		f4 = mtrx[re][5];
		if (f4=="I" || f4=="O")
		{
			f4="S";
		}
		
		f=mtrx[re][3] + mtrx[re][1] + f4 + "HH";
		
        A.ans1 = f;
        A.ans2 = mtrx[re][7];
        A.ans3 = mtrx[re][9];
        A.symb1 = "" ;
        A.symb2 = "" ; 
        /* assemble the answer for this condition */
		

		A=sul(A,2); // r=2 to give correct Bessel function

		Symb1 += A.symb1;
		Symb2 += A.symb2;
		Eigans1 += A.ans1;
		Eigans2 += "\\\\ \n" + A.ans2;
		Eigans53 += "\\;" + A.ans5 + "\\;" + A.ans3;
		Eigans4 += A.ans4 + "\\\\";
	// That builds the answer string for this factor.
	}
// Build nonhomf
A = KAPPAGR(A);

//   Write to outfile
outFilE << Symb1 << Eigans1 << "\\;" << A.AAA << Eigans2  << "\\\\" ;
outFilE << A.AAA << Symb2 << Eigans53 << "\\;" << A.nonhomf << Eigans4 ;
 return A;
}

//  *********  KandF  ********


pdeans KandF(pdeans C, std::string BC, int r)
//  Puts the nonhomogeneous parameter string in KAPPA, the coeff in AAA, and the nonhomog function in nonhomf
{
std::string KAPPA, nonhomf, AAA;

	if(eqtype[0]==0)   //  LaPlace
	{
		if(eqtype[1]==2)  //  2-D
		{
			if(r==0)  // x
			{
				KAPPA = "\\kappa_y";
				AAA="A(\\kappa_y)";
				nonhomf="f_{x=" + BC + "}(y)";
			}
			if(r==1)  // y
			{
				KAPPA = "\\kappa_x";
				AAA="A(\\kappa_x)";
				nonhomf="f_{y=" + BC + "}(x)";
			}
			if(r==3)   // theta
			{
				KAPPA = "\\kappa_r";
				AAA="A(\\kappa_r)";
				nonhomf="f_{\\theta=" + BC + "}(r)";
			}
			if(r==4) //   r 2_D
			{
				KAPPA = "\\kappa_{\\theta}";
				AAA="A(\\kappa_{\\theta})";
				nonhomf="f_{r=" + BC + "}(\\theta)";
			}
		}
		if(eqtype[1]==3)   //  3_D
		{
			if(r==0)  // x
			{
				KAPPA="\\sqrt{\\kappa_y^2+\\kappa_z^2}";
				AAA="A(\\kappa_y, \\, \\kappa_z)";
				nonhomf="f_{x=" + BC + "}(y,z)";
			}
			if(r==1)  //y
			{
				KAPPA="\\sqrt{\\kappa_x^2+\\kappa_z^2}";
				AAA="A(\\kappa_x, \\, \\kappa_z)";
				nonhomf="f_{y=" + BC + "}(x,z)";
			}
			if(r==2)  //z
			{
				if(eqtype[2]==0)  // rect
				{
					KAPPA="\\sqrt{\\kappa_x^2+\\kappa_y^2}";
					AAA="A(\\kappa_x, \\, \\kappa_y)";
					nonhomf="f_{z=" + BC + "}(x,y)";
				}
				if(eqtype[2]==1)  // cyl
				{
					KAPPA="\\kappa_{r:\\theta}";
					AAA = "A(\\kappa_{\\theta},\\kappa_{r:\\theta})";
					nonhomf="f_{z=" + BC + "}(\\theta,r)";
				}
			}

			if(r==3)   // theta
			{
				KAPPA="\\kappa_{r:z}";
				AAA = "A(\\kappa_z,\\kappa_{r:z})";
				nonhomf="f_{\\theta=" + BC + "}(z,r)";
			}
			if(r==5)    //  rho
			{
				KAPPA = "\\kappa_{\\theta},\\kappa_z";  //  not used
				AAA = "A(\\kappa_{\\theta},\\kappa_z)";
				nonhomf="f_{r=" + BC + "}(\\theta, z)";
			}
			if(r==6)   // r 3-D
			{
				KAPPA = "\\ell m";   // not used
				AAA="A_{\\ell m}";
				nonhomf="f_{r=" + BC + "}(\\theta,\\phi)";
			}
		}
	}

	if(eqtype[0]==1)   //  diffusion, time
	{
		if(eqtype[1]==2)  // 2_D  (1-D is special)
			{
			if(r==0)  //x
			{
				KAPPA = "\\kappa_y";
				AAA="A(\\kappa_y)";
				nonhomf="f_{x=" + BC + "}(y)";
			}
			if(r==1)  //y
			{
				KAPPA = "\\kappa_x";
				AAA="A(\\kappa_x)";
				nonhomf="f_{y=" + BC + "}(x)";
			}
			if(r==3)  // theta
			{
				KAPPA = "\\kappa_{r}";
				AAA="A(\\kappa_{r})";
				nonhomf="f_{\\theta=" + BC + "}(r)";
			}
			if(r==4)  // r 2-D   
			{
				KAPPA = "\\kappa_{\\theta}";
				AAA="A(\\kappa_{\\theta})";
				nonhomf="f_{r=" + BC + "}(\\theta)";
			}
		}
		if(eqtype[1]==3)   //  3-D
		{
			if(r==0)    // x
			{
				KAPPA="\\sqrt{\\kappa_y^2+\\kappa_z^2}";
				AAA="A(\\kappa_y, \\, \\kappa_z)";
				nonhomf="f_{x=" + BC + "}(y,z)";
			}
			if(r==1)  // y
			{
				KAPPA="\\sqrt{\\kappa_x^2+\\kappa_z^2}";
				AAA="A(\\kappa_x, \\, \\kappa_z)";
				nonhomf="f_{y=" + BC + "}(x,z)";
			}
			if(r==2)  // z
			{
				if(eqtype[2]==0)  // rect
				{
					KAPPA="\\sqrt{\\kappa_x^2+\\kappa_y^2}";
					AAA="A(\\kappa_x, \\, \\kappa_y)";
					nonhomf="f_{z=" + BC + "}(x,y)";
				}
				if(eqtype[2]==1)   // cyl
				{
					KAPPA="\\kappa_{r:\\theta}";
					AAA = "A(\\kappa_{\\theta},\\kappa_{r:\\theta})";
					nonhomf="f_{z=" + BC + "}(\\theta,r)";
				}
			}

			if(r==3)  // theta
			{
				KAPPA="\\kappa_{r:z}";
				AAA = "A(\\kappa_z,\\kappa_{r:z})";
				nonhomf="f_{\\theta=" + BC + "}(z,r)";
			}
			if(r==5)   //  rho
			{
				KAPPA="\\kappa_z";
				AAA = "A(\\kappa_{\\theta},\\kappa_z)";
				nonhomf="f_{r=" + BC + "}(\\theta, z)";
			}
			if(r==6)   //  r 3-D
			{
				KAPPA = "\\ell m";   // not used
				AAA="A_{\\ell m}";
				nonhomf="f_{r=" + BC + "}(\\theta,\\phi)";
			}
		}
	}
	if(eqtype[0]==2)    //diffusion, s plane
	{
		if(eqtype[1]==1)   // 1-D
			{
				KAPPA="\\sqrt{s}";
				AAA="A(s)";
				nonhomf="F_{x=" + BC + "}(s)";
			}
		
		
		if(eqtype[1]==2)   // 2-D
		{
			if(r==0)   // x
			{
				KAPPA="\\sqrt{\\kappa_y^2+s}";
				AAA="A(s;\\kappa_y)";
				nonhomf="F_{x=" + BC + "}(s;y)";
			}
			if(r==1)   //y
			{
				KAPPA="\\sqrt{\\kappa_x^2+s}";
				AAA="A(s;\\kappa_x)";
				nonhomf="F_{y=" + BC + "}(s;x)";
			}
			if(r==3)  // theta   
			{
				KAPPA = "\\sqrt{\\kappa_r^2+s}";
				AAA="A(s;\\kappa_r)";
				nonhomf="f_{\\theta=" + BC + "}(s;r)";
			}
			if(r==4)  // r 2-D   
			{
				KAPPA = "\\sqrt{s}";
				AAA="A(s;\\kappa_{\\theta})";
				nonhomf="f_{r=" + BC + "}(s;\\theta)";
			}  //**************???????????????
		}
		if(eqtype[1]==3)   // 3-D
		{
			if(r==0)   // x
			{
				KAPPA="\\sqrt{\\kappa_y^2+\\kappa_z^2+s}";
				AAA="A(s;\\kappa_y,\\kappa_z)";
				nonhomf="F_{x=" + BC + "}(s;y,z)";
			}
			if(r==1)   // y
			{
				KAPPA="\\sqrt{\\kappa_x^2+\\kappa_z^2+s}";
				AAA="A(s;\\kappa_x,\\kappa_z)";
				nonhomf="F_{y=" + BC + "}(s;x,z)";
			}
			if(r==2)   // z
			{
			if(eqtype[2]==0)   // rect
				{
				KAPPA="\\sqrt{\\kappa_x^2+\\kappa_y^2+s}";
				AAA="A(s;\\kappa_x,\\kappa_y)";
				nonhomf="F_{z=" + BC + "}(s;x,y)";
				}
			if(eqtype[2]==1)   // cyl
				{
				KAPPA="\\sqrt{\\kappa_{r : \\theta}^2+s}";