#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#define e 0.0018
#define t 3500.  

double gauss(double (*f)(),double a,double b,int n)
{
	double c,d,s;
	double	y[5],h[5];
	int i;
	c=(b-a)/2;
	d=(b+a)/2;
	switch(n)
	{
	case(1):
		y[0]=0;  	h[0]=2;
		break;
	case(2): 
		y[0]=-1/sqrt(3);    h[0]=1;
		y[1]=1/sqrt(3); 	h[1]=1;
		break;
	case(3):
		y[0]=-sqrt(15)/5;   h[0]=5.0/9;
		y[1]=0;             h[1]=8.0/9;
		y[2]=sqrt(15)/5;	h[2]=5.0/9;
		break;
	case(4):
		y[0]=0.8611363;   h[0]=0.3478548;
		y[1]=0.3398810;	  h[1]=0.6521452;
		y[2]=-0.8611363;  h[2]=0.3478548;
		y[3]=-0.3398810;  h[3]=0.6521452;
		break;
	default:
		y[0]=0.9061793;	  h[0]=0.2369269;
		y[1]=0.5384693;   h[1]=0.4786287;
		y[2]=0;           h[2]=0.5688889;
		y[3]=-0.9061793;  h[3]=0.2369269;
		y[4]=-0.5384693;  h[4]=0.4786287;
    	break;
	}
	for(i=0,s=0;  i<n;  i++)
		s+=h[i]*f(c*y[i]+d,b);
	s*=c;
	return(s);
}

double f(double y, double b)
{
 	return((0.131*e*y/b-20*e*e*y*y/(b*b))*400/(0.00000484+0.002155*e*y/b));  //計算 積分σb dy
 }
double f1(double y, double b)
{
 	return((0.131*e*y/b-20*e*e*y*y/(b*b))*400*(560-b+y)/(0.00000484+0.002155*e*y/b)); //計算積分σb（h0-x+y) dy
 }

void main(void)
{
	FILE *fp,*fo;
	double a,b,b0,y,N,M,z,s1,s2,s;   //其中s1為受拉鋼筋的合力;s2為受壓鋼筋的合力;s為鋼筋屈服時的合力
	double	*x,*num;
	double b1,*xh,NO[2];
	int n,i,j,k;
	a=0.0;
	n=5;
	s=942*455;
    xh=(double*)malloc(601*sizeof(double));
	for(i=0;i<=600;i++) 
		xh[i]=i;	
	for(j=0; j<600;j++)
	{
		
		for(k=0; k<2; k++)
		{
			b1=xh[j];
			y=gauss(f,a,b1,n);
			z=gauss(f1,a,b1,n);
			s1=200000*e*(560-b1)*942.48/b1;
			s2=200000*e*(b1-40)*942.48/b1;
			if(s2>=s) s2=s;
			if(s1>=s) s1=s;
			NO[k]=(y+s2-s1)/1000;          //計算軸力N
			j++;
		}
		j=j-2;
		if(NO[0]<=t && NO[1]>t)
		{
			b0=xh[j];
			break;
		}
	}
	fp=fopen("datas.text","w");
	num=(double*)malloc(40001*sizeof(double));
	num[0]=b0;
	for(k=0;k<=40000;k++)
	{
		num[k+1]=num[k]+0.000025;
		fprintf(fp,"%.10lf\n",num[k]);
	}
	fclose(fp);
	fo=fopen("datas.text","r");
	x=(double*)malloc(40001*sizeof(double));
	for(i=0; i<=40000; i++ )
		fscanf(fo,"%lf",x+i);
	fclose(fo);
	for(i=0; i<=40000;i++)
	{
		b=x[i];
		y=gauss(f,a,b,n);   
		z=gauss(f1,a,b,n);   
		s2=200000*e*(b-40)*942.48/b;
		s1=200000*e*(560-b)*942.48/b;
		if(s2>=s) s2=s;
		if(s1>=s) s1=s;
 	    N=(y+s2-s1)/1000;      //計算軸力N
     	if(fabs(N-t)<=0.0002) 
		{
				z=gauss(f1,a,b,n);
                M=(z+s2*(560-40)-t*1000*260)/1000000; //計算彎矩Ne
				printf("e=%lf\n\nb=%lf\nN=%lf\nM=%lf\ns1=%lf\ns2=%lf\n",e,b,N,M,s1/942.48,s2/942.48);
		    break;
		}
		
	}
    	
		
}