import java.util.*;
public class HT_infinite_wall2
{
	//incomplete gamma function P(a,x)
public static double gammp(double a,double x)
{
	double gamser,gammcf;
	if (x < 0.0 || a <= 0.0) System.out.println("Invalid arguments in routine GAMMP");
	if (x < (a+1.0)) {
		gamser=gser(a,x);
		return gamser;
	} else {
		gammcf=gcf(a,x);
		return 1.0-gammcf;
	}
}

//incomplete gamma function Q(a,x)=1-P(a,x);
public static double gammq(double a,double x)
{
	double gamser,gammcf;
	if (x < 0.0 || a <= 0.0) System.out.println("Invalid arguments in routine GAMMQ");
	if (x < (a+1.0)) {
		gamser=gser(a,x);
		return 1.0-gamser;
	} else {
		gammcf=gcf(a,x);
		return gammcf;
	}
}

public static double gser(double a,double x)
{
	int ITMAX=200;
	double EPS=3.0e-7;
	int n;
	double sum,del,ap;
	double gln=lngamma(a)[0];
	double gamser=0;
	if (x <= 0.0) {
		if (x < 0.0) System.out.println("x less than 0 in routine GSER");
		gamser=0.0;
		return gamser;
	} else {
		ap=a;
		del=sum=1.0/a;
		for (n=1;n<=ITMAX;n++) {
			ap += 1.0;
			del *= x/ap;
			sum += del;
			if (Math.abs(del) < Math.abs(sum)*EPS) {
				gamser=sum*Math.exp(-x+a*Math.log(x)-(gln));
				return gamser;
			}
		}
		//System.out.println("a too large, ITMAX too small in routine GSER");
		return gamser;
	}
}

public static double gcf(double a,double x)
{
	int ITMAX=200;
	double EPS=3.0e-7;
	int n;
	double gln;
	double gammcf=0;
	double gold=0.0,g,fac=1.0,b1=1.0;
	double b0=0.0,anf,ana,an,a1,a0=1.0;
	gln=lngamma(a)[0];
	a1=x;
	for (n=1;n<=ITMAX;n++) {
		an=(double) n;
		ana=an-a;
		a0=(a1+a0*ana)*fac;
		b0=(b1+b0*ana)*fac;
		anf=an*fac;
		a1=x*a0+anf*a1;
		b1=x*b0+anf*b1;
		if (a1!=0) {
			fac=1.0/a1;
			g=b1*fac;
			if (Math.abs((g-gold)/g) < EPS) {
				gammcf=Math.exp(-x+a*Math.log(x)-gln)*g;
				return gammcf;
			}
			gold=g;
		}
	}
	//System.out.println("a too large, ITMAX too small in routine GCF");
    return gammcf;
}
//error function
public static double  erf(double x)
{
   //error function
   return x < 0.0 ? -gammp(0.5,x*x) : gammp(0.5,x*x);
}
//natural logarithm of gamma function 
	public static double [] lngamma(double x)
    {//logaithm of gamma function
    //xx[0] value of logarithm of gamma function 
    //xx[1] sign of gamma functon
	double xx[]=new double[2];    
    double result;
    double a;
    double b;
    double c;
    double p;
    double q;
    double u;
    double w;
    double z;
    int i;
    double logpi;
    double ls2pi;
    double tmp;
    double sgngam = 1;
    logpi = 1.14472988584940017414;
    ls2pi = 0.91893853320467274178;
    if( x<-34.0 )
    {
        q = -x;
        double yy[]=lngamma(q);
        w = yy[0];
        tmp=yy[1];
        p = Math.floor(q);
        i = (int)Math.round(p);
        if( i%2==0 )
        {
            sgngam = -1;
        }
        else
        {
            sgngam = 1;
        }
        z = q-p;
        if( z>0.5 )
        {
            p = p+1;
            z = p-q;
        }
        z = q*Math.sin(Math.PI*z);
        result = logpi-Math.log(z)-w;
        xx[0]=result;
        xx[1]=sgngam;
        return xx;
    }
    if( x<13 )
    {
        z = 1;
        p = 0;
        u = x;
        while(u>=3)
        {
            p = p-1;
            u = x+p;
            z = z*u;
        }
        while(u<2)
        {
            z = z/u;
            p = p+1;
            u = x+p;
        }
        if( z<0 )
        {
            sgngam = -1;
            z = -z;
        }
        else
        {
            sgngam = 1;
        }
        if( u==2 )
        {
            result = Math.log(z);
            //return result;
            xx[0]=result;
            xx[1]=sgngam;
            return xx;
        }
        p = p-2;
        x = x+p;
        b = -1378.25152569120859100;
        b = -38801.6315134637840924+x*b;
        b = -331612.992738871184744+x*b;
        b = -1162370.97492762307383+x*b;
        b = -1721737.00820839662146+x*b;
        b = -853555.664245765465627+x*b;
        c = 1;
        c = -351.815701436523470549+x*c;
        c = -17064.2106651881159223+x*c;
        c = -220528.590553854454839+x*c;
        c = -1139334.44367982507207+x*c;
        c = -2532523.07177582951285+x*c;
        c = -2018891.41433532773231+x*c;
        p = x*b/c;
        result = Math.log(z)+p;
        //return result;
        xx[0]=result;
        xx[1]=sgngam;
        return xx;
    }
    q = (x-0.5)*Math.log(x)-x+ls2pi;
    if( x>100000000 )
    {
        result = q;
        //return result;
        xx[0]=result;
        xx[1]=sgngam;
        return xx;
    }
    p = 1/(x*x);
    if( x>=1000.0 )
    {
        q = q+((7.9365079365079365079365*0.0001*p-2.7777777777777777777778*0.001)*p+0.0833333333333333333333)/x;
    }
    else
    {
        a = 8.11614167470508450300*0.0001;
        a = -5.95061904284301438324*0.0001+p*a;
        a = 7.93650340457716943945*0.0001+p*a;
        a = -2.77777777730099687205*0.001+p*a;
        a = 8.33333333333331927722*0.01+p*a;
        q = q+a/x;
    }
    xx[0]= q;
    xx[1]=sgngam;
    return xx;
    }
//complimentary error function
public static double  erfc(double x)
{
	return x < 0.0 ? 1.0+gammp(0.5,x*x) : gammq(0.5,x*x);
}
public static double T(double x,double t,double Tinf,double alpha,double H)
{
//H = h/k
double x1=x/Math.sqrt(4.0*alpha*t);
double TT=Tinf*(erfc(x1)-Math.exp(H*x+H*H*alpha*t)*erfc(H*Math.sqrt(alpha*t)+x1));
return TT; 	
}
     public static void main (String args[])  
     {  
       double alpha=4.0e-5;
       double H=2.0;
       double x=0.1; //m
       double t=300.0;//s
       double T1=20.0;// degree C;
       double Tinf=100.0;
       double T=T(x,t,Tinf,alpha,H)+T1;
       System.out.println("T="+T);
     }
}