/* composite.c calulate composite expansion for 2-3 or 2-2 case
see bleow about commenting line to reduce 2-3 to 2-2 expansion
*/
#include <stdio.h>
#include <math.h>
double h0[10001],dh0[10001],ddh0[10001],h,a0,b0,df,e;
double h1[10001],dh1[10001],ddh1[10001],a2,b2,c2,eta,g2inf;
double g0[10001],dg0[10001],ddg0[10001],alpha0,sqrt_alpha0;
double fa,fb,fc,fd,dfa,dfb,dfc,dfd,ddfa,ddfb,ddfc,ddfd;
double ga,gb,gc,gd,dga,dgb,dgc,dgd,ddga,ddgb,ddgc,ddgd;
double sp,sm,mu0,mu1,gamma0,gamma1,beta0;
double lambda0,lambda1;
double u[101][21],x,y,dx,dy,ny,nx,s,eta,sqrt2,xthird,dddfa;
double xpts[21];
int i,j,k,l,m,n,k1,k2;
FILE *fp;
main()
{
	sqrt2=pow((double)2.,(double)0.5);

/* calculate base flow to eta=10, s=10 */

		
		h=0.01;
		n=2000;
		g0[0]  =0.;
		dg0[0] =0.;
		ddg0[0]=1.;
		for(i=1;i<=n;i++)
			{
			fa=   h*dg0[i-1];
			dfa=  h*ddg0[i-1];
			ddfa=-h*(g0[i-1]*ddg0[i-1]);

			fb=   h*(dg0[i-1]+dfa/2);
			dfb=  h*(ddg0[i-1]+ddfa/2);
			ddfb=-h*((g0[i-1]+fa/2)*(ddg0[i-1]+ddfa/2));

			fc=   h*(dg0[i-1]+dfb/2);
			dfc=  h*(ddg0[i-1]+ddfb/2);
			ddfc=-h*((g0[i-1]+fb/2)*(ddg0[i-1]+ddfb/2)) ;

			fd=   h*(dg0[i-1]+dfc);
			dfd=  h*(ddg0[i-1]+ddfc);
			ddfd=-h*((g0[i-1]+fc)*(ddg0[i-1]+ddfc));


			g0[i]=g0[i-1]+(fa+2*fb+2*fc+fd)/6;
			dg0[i]=dg0[i-1]+(dfa+2*dfb+2*dfc+dfd)/6;
			ddg0[i]=ddg0[i-1]+(ddfa+2*ddfb+2*ddfc+ddfd)/6;
			}
		alpha0=dg0[n];
		sqrt_alpha0=pow(alpha0,(double)0.5);
		beta0=0.;
		for(i=1;i<=n;i++)
			{
			beta0=beta0+0.5*h*(2*alpha0-dg0[i]-dg0[i-1]);
			}
		beta0=beta0/sqrt_alpha0;
		fp=fopen("bl.out","w");
		fprintf(fp,"data\n");
		for(i=0;i*h*sqrt_alpha0<=5.;i+=10)
			{
			fprintf(fp,"%12.4lf %12.5le\n",i*h*sqrt_alpha0,dg0[i]/alpha0);
			}
		fclose(fp);

		h0[0]  =0.;
		dh0[0] =1.;
		ddh0[0]=0.;
		for(i=1;i<=n;i++)
			{
			fa=   h*dh0[i-1];
			dfa=  h*ddh0[i-1];
			ddfa= h*(dh0[i-1]*dh0[i-1]-2*h0[i-1]*ddh0[i-1])/3;

			fb=   h*(dh0[i-1]+dfa/2);
			dfb=  h*(ddh0[i-1]+ddfa/2);
			ddfb= h*((dh0[i-1]+dfa/2)*(dh0[i-1]+dfa/2)-
                                   2*(h0[i-1]+fa/2)*(ddh0[i-1]+ddfa/2))/3;

			fc=   h*(dh0[i-1]+dfb/2);
			dfc=  h*(ddh0[i-1]+ddfb/2);
			ddfc= h*((dh0[i-1]+dfb/2)*(dh0[i-1]+dfb/2)-
                                   2*(h0[i-1]+fb/2)*(ddh0[i-1]+ddfb/2))/3;

			fd=   h*(dh0[i-1]+dfc);
			dfd=  h*(ddh0[i-1]+ddfc);
			ddfd= h*((dh0[i-1]+dfc)*(dh0[i-1]+dfc)-
                                   2*(h0[i-1]+fc)*(ddh0[i-1]+ddfc))/3;


			h0[i]=h0[i-1]+(fa+2*fb+2*fc+fd)/6;
			dh0[i]=dh0[i-1]+(dfa+2*dfb+2*dfc+dfd)/6;
			ddh0[i]=ddh0[i-1]+(ddfa+2*ddfb+2*ddfc+ddfd)/6;
			}
		for(i=n-20;i<=n;i++)
			{				
			fc=i*h;
			gamma0=ddh0[i]/2;
			fa=h0[i]/gamma0;
			fb=pow(fa,(double) 0.5);
			mu0=fb-fc;

			}

		h1[0]=0.;
		dh1[0]=1.;
		ddh1[0]=0.;
		for(i=1;i<=n;i++)
			{
			fa=(h0[i-1]+h0[i])/2;
			dfa=(dh0[i-1]+dh0[i])/2;
			ddfa=(ddh0[i-1]+ddh0[i])/2;

			ga=   h*dh1[i-1];
			dga=  h*ddh1[i-1];
			ddga= h*(5*dh0[i-1]*dh1[i-1]-5*ddh0[i]*h1[i-1]-2*h0[i-1]*ddh1[i-1])/3;

			gb=   h*(dh1[i-1]+dga/2);
			dgb=  h*(ddh1[i-1]+ddga/2);
			ddgb= h*(5*dfa*(dh1[i-1]+dga/2)-5*ddfa*(h1[i-1]+ga/2)-
            					2*fa*(ddh1[i-1]+ddga/2))/3;
			gc=   h*(dh1[i-1]+dgb/2);
			dgc=  h*(ddh1[i-1]+ddgb/2);
			ddgc= h*(5*dfa*(dh1[i-1]+dgb/2)-5*ddfa*(h1[i-1]+gb/2)-
            					2*fa*(ddh1[i-1]+ddgb/2))/3;

			gd=   h*(dh1[i-1]+dgc);
			dgd=  h*(ddh1[i-1]+ddgc);
			ddgd= h*(5*dh0[i]*(dh1[i-1]+dgc)-5*ddh0[i]*(h1[i-1]+gc)-
            					2*h0[i]*(ddh1[i-1]+ddgc))/3;
			h1[i]=h1[i-1]+(ga+2*gb+2*gc+gd)/6;
			dh1[i]=dh1[i-1]+(dga+2*dgb+2*dgc+dgd)/6;
			ddh1[i]=ddh1[i-1]+(ddga+2*ddgb+2*ddgc+ddgd)/6;
			}
		for(i=n-20;i<=n;i++)
			{
			sp=i*h+mu0;
			dfa=sp*sp*(sp*sp*sp+30/gamma0);
			gamma1=ddh1[i]/(20*sp*sp*sp+60./gamma0);
			mu1=(h1[i]/gamma1-dfa)/sp;				
			
			}
		fa=4.*gamma0*pow(alpha0,(double)1.5)/sqrt2;
		fa=1./fa;
		lambda0=pow(fa,(double)(1./3.));
		fa=pow((double)2.,(double)2.5)*120.*gamma1*pow(alpha0,(double)3.)*
							pow(lambda0,(double)5.);
		lambda1=-1./fa;

	printf("Constant list:\n");
	printf("beta_0   %15.8le\n",beta0);
	printf("alpha_0  %15.8le\n",alpha0);
	printf("gamma_0  %15.8le\n",gamma0);
	printf("mu_0     %15.8le\n",mu0);
	printf("lambda_0 %15.8le\n",lambda0);
	printf("gamma_1  %15.8le\n",gamma1);
	printf("mu_1     %15.8le\n",mu1);
	printf("lambda_1 %15.8le\n",lambda1);
	xpts[0]=0.;
	xpts[1]=.01;
	xpts[2]=0.05;
	xpts[3]=0.1;
	xpts[4]=0.2;
	xpts[5]=0.3;
	xpts[6]=0.4;
	xpts[7]=0.5;
	xpts[8]=0.5;
	xpts[9]=0.5;
	xpts[10]=0.75;
	xpts[11]=1.;
	dx=0.01;
	dy=.05;
	nx=7;
	ny=100;
	i=0;
	for(j=0;j<=ny;j++)
		{
		y=j*dy;eta=y/(sqrt2*sqrt_alpha0);k1=eta/h;fa=(eta-k1*h)/h;
		u[j][0]=((1-fa)*dg0[k1]+fa*dg0[k1+1])/alpha0;
		}
	for(i=1;i<=nx;i++)
		{
		for(j=0;j<=ny;j++)
			{
			x=xpts[i];
			xthird=pow(x,(double)(1./3.));
			y=j*dy;
			s=lambda0*y/xthird;
			eta=y/(sqrt2*sqrt_alpha0);
			k1=eta/h;
			fa=(eta-k1*h)/h;
			k2=s/h;
			fb=(s-k2*h)/h;
			dfa=((1-fa)*dg0[k1]+fa*dg0[k1+1])/alpha0;
			ddfa=((1-fa)*ddg0[k1]+fa*ddg0[k1+1])/(alpha0*sqrt_alpha0);
			dddfa=-((1-fa)*g0[k1]+fa*g0[k1+1])*((1-fa)*ddg0[k1]+
				fa*ddg0[k1+1])/(alpha0*alpha0);
			dga=lambda0*lambda0*((1-fb)*dh0[k2]+fb*dh0[k2+1]);
			dgb=lambda0*lambda1*((1-fb)*dh1[k2]+fb*dh1[k2+1]);
			s=lambda0*y+mu0*xthird;
			u[j][i]=
				xthird*(dga+x*dgb)
				+dfa+xthird*mu0*ddfa/(sqrt2*lambda0)
/* comment out the next line to get 2-2 xomposite output */
		  +mu0*mu0*xthird*xthird*dddfa/(4*sqrt2*lambda0*lambda0)
				-2*lambda0*lambda0*gamma0*s
		-lambda0*lambda1*gamma1*(5*s*s*s*s+60.*x*s/gamma0+x*xthird*mu1) 
				;
			}
		}

	
		fp=fopen("composite.out","w");
		fprintf(fp,"y max:1.0\n");
		for(i=0;i<=nx;i++)
			{
			fprintf(fp,"data\n");
			for(j=0;j<=ny;j++)
			{fprintf(fp,"%8.2lf %12.6le\n",j*dy,u[j][i]);
			}
			}
		fclose(fp);
				
}