/*
calculate limiting values of f_2 and c_{4,2} to show failure
of simple power series expansion
*/
#include <stdio.h>
#include <math.h>
double f0[10001],df0[10001],ddf0[10001],h,a0,b0,df,e,r,da2,a2offset;
double f2[10001],df2[10001],ddf2[10001],a2,b2,c2,eta,g2inf;
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 alpha0,sqrt_alpha0,a3t;
double f2lim[100],c2lim[100];
int i,j,k,l,m,n,kpoints;
FILE *fp;
main()
{
/* calculate base flow */
	alpha0=(double)1.65519036;
	sqrt_alpha0=pow(alpha0,(double)0.5);
	a3t=pow(alpha0,(double)(-1.5));
	h=0.002;
	n=5000;
	f0[0]=0.;
	df0[0]=0.;
	ddf0[0]=a3t;
	for(i=1;i<=n;i++)
		{
		fa=   h*df0[i-1];
		dfa=  h*ddf0[i-1];
		ddfa= -h*f0[i-1]*ddf0[i-1];
		fb=   h*(df0[i-1]+dfa/2);
		dfb=  h*(ddf0[i-1]+ddfa/2);
		ddfb= -h*(f0[i-1]+fa/2)*(ddf0[i-1]+ddfa/2);
		fc=   h*(df0[i-1]+dfb/2);
		dfc=  h*(ddf0[i-1]+ddfb/2);
		ddfc= -h*(f0[i-1]+fb/2)*(ddf0[i-1]+ddfb/2);
		fd=   h*(df0[i-1]+dfc);
		dfd=  h*(ddf0[i-1]+ddfc);
		ddfd= -h*(f0[i-1]+fc)*(ddf0[i-1]+ddfc);

		f0[i]=f0[i-1]+(fa+2*fb+2*fc+fd)/6;
		df0[i]=df0[i-1]+(dfa+2*dfb+2*dfc+dfd)/6;
		ddf0[i]=ddf0[i-1]+(ddfa+2*ddfb+2*ddfc+ddfd)/6;
		}
	b0=0.;
	for(i=1;i<=n;i++)
		{
		b0=b0+0.5*h*(2*df0[n]-df0[i-1]-df0[i]);
		}
	printf("b0= %12.5le\n",b0);
	printf("input da2, kpoints,a2offset\n");
	scanf("%lf",&da2);
	scanf("%d",&kpoints);
	scanf("%lf",&a2offset);
	for(k=0;k<=kpoints;k++)
		{
		a2=k*da2-a2offset;
		f2[0]=0.;
		df2[0]=0.;
		ddf2[0]=a2*a3t;
		for(i=1;i<=n;i++)
			{
			e=(i-1)*h;
			r=(f0[i-1]-e*df0[i-1])*(f0[i-1]-e*df0[i-1])-b0*b0;
			fa=   h*df2[i-1];
			dfa=  h*ddf2[i-1];
			ddfa= h*(r-f0[i-1]*ddf2[i-1]
			  	-2*df0[i-1]*df2[i-1]
			 	 +ddf0[i-1]*f2[i-1]);

			e=e+0.5*h;
			ga=0.5*(f0[i-1]+f0[i]);
			dga=0.5*(df0[i-1]+df0[i]);
			ddga=0.5*(ddf0[i-1]+ddf0[i]);
			r=(ga-e*dga)*(ga-e*dga)-b0*b0;

			fb=   h*(df0[i-1]+dfa/2);
			dfb=  h*(ddf0[i-1]+ddfa/2);
			ddfb= h*(r-ga*(ddf2[i-1]+ddfa/2)
				 -2*dga*(df2[i-1]+dfa/2)
				 +ddga*(f2[i-1]+fa/2));

			fc=   h*(df0[i-1]+dfb/2);
			dfc=  h*(ddf0[i-1]+ddfb/2);
			ddfc= h*(r-ga*(ddf2[i-1]+ddfb/2)
				 -2*dga*(df2[i-1]+dfb/2)
				 +ddga*(f2[i-1]+fb/2));

			e=e+0.5*h;
			r=(f0[i]-e*df0[i])*(f0[i]-e*df0[i])-b0*b0;
			fd=   h*(df0[i-1]+dfc);
			dfd=  h*(ddf0[i-1]+ddfc);
			ddfd= h*(r-f0[i]*(ddf2[i-1]+ddfc)
			 	 -2*df0[i]*(df2[i-1]+dfc)
			 	 +ddf0[i]*(f2[i-1]+fc));

			f2[i]=f2[i-1]+(fa+2*fb+2*fc+fd)/6;
			df2[i]=df2[i-1]+(dfa+2*dfb+2*dfc+dfd)/6;
			ddf2[i]=ddf2[i-1]+(ddfa+2*ddfb+2*ddfc+ddfd)/6;
			}
		c2=0.;
		for(i=1;i<=n;i++)
			{
			eta=i*h;
			c2=c2+h*(eta*df2[i]+f2[i]-f2[n]);
			}
		f2lim[k]=f2[n];
		c2lim[k]=c2;
			
		}
	fp=fopen("test.out","w");
	fprintf(fp,"data\n");
	for(i=0;i<=kpoints;i++)
		{
		fprintf(fp,"%8.4lf %15.7le\n",i*da2-a2offset,f2lim[i]);
		}
	fprintf(fp,"data\n");
	for(i=0;i<=kpoints;i++)
		{
		fprintf(fp,"%8.4lf %15.7le\n",i*da2-a2offset,c2lim[i]);
		}
	fclose(fp);
}