#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#define NR 8
#define NC 9

static double A[NR][NC], a[NR+1], a0[NR+1];

double g(double x, double a, double b)		/* x - phi, a - g, b - theta */
{
	double z;
	z=a*sin(2*(b+x))-sin(x);
	return z;
}


double g2(double x, double a, double b)
{
	double z;
	z=-4*a*sin(2*(b+x))+sin(x);
	return z;
}


double aniso_magn(double a, double b)		/* x - phi, a - g, b - theta */
{
	int i, j, n;
	double dx, f1, f, x, x1, y, z;
	double alpha, beta, pi;
	
	pi=M_PI;
	n=50;
	dx=pi/4/n;
	while(b>pi/2) b-=pi;
	while(b<-pi/2) b+=pi;
	if(b<0)
	{
		x=0;
		while(g(x, a, b)<0)
		{
			x-=dx;
		}
		alpha=x;
		while(g2(x, a, b)>0)
		{
			x-=dx;
		}
		beta=x;
	}
	else
	{
		x=0;
		while(g(x, a, b)>0)		/* x - phi, a - g, b - theta */
		{
			x+=dx;
		}
		alpha=x;
		while(g2(x, a, b)<0)		/* second derivative of g(x) */
		{
			x+=dx;
		}
		beta=x;
	}
	x=(alpha+beta)/2;
	z=1.0;
	j=0;
	while(fabs(z)>0.001 && j<10)
	{
		f=a*sin(2*(x+b))-sin(x);
		f1=2*a*cos(2*(x+b))-cos(x);
		x1=x-f/f1;
		if(f1==0) {x=0; break;}
		z=x1-x;
		x=x1;
		j++;
	}
	return x;
}


double aniso_magn_r(double a, double b)		/* x - phi, a - g, b - theta */
{
	char dummy[30];
	int i, j, n;
	double dx, f1, f, x, x1, y, z;
	double pi;
	
	pi=M_PI;
	n=50;
	dx=pi/4/n;
	if(a<0){a*=-1; b+=pi/2;}
	while(b>=pi/2) b-=pi;
	while(b<-pi/2) b+=pi;
	if(a>0.5 && b<0)
	{ 
		x=-pi/2;
		while(g2(x,a,b)<0) x+=dx;
	}
	else if(a>0.5 && b>=0)
	{
		x=pi/2;
		while(g2(x,a,b)>0) x-=dx;
	}
	else if(b<0)
	{
		x=-pi/2;
		while(g2(x, a, b)<0)
		{
			x+=dx;
		}
	}
	else
	{
		x=pi/2;
		while(g2(x, a, b)>0)
		{
			x-=dx;
		}
	}
	z=1.0;
	j=0;
	while(fabs(z)>0.001 && j<10)
	{
		f=a*sin(2*(x+b))-sin(x);
		f1=2*a*cos(2*(x+b))-cos(x);
		x1=x-f/f1;
		if(f1==0) {x=0; break;}
		z=x1-x;
		x=x1;
		j++;
	}
	return x;
}


int LeastSquare(int n)
{
	char buffer[10];
	long double  B[NR], pivot, Z, biggest;
	int i, j, k, ipivot, nr, nc;
	
	nr=n;
	nc=nr+1;
	
	for(j=0;j<nr;j++){
		biggest=-1;
		for(i=j;i<nr;i++){
//			printf("%lf\n", A[i][j]);
			if(fabs(A[i][j])>biggest){
				biggest=fabs(A[i][j]);
				ipivot=i;
			}
		}
		if(ipivot!=j){
			for(k=0;k<nc;k++){
				B[k]=A[ipivot][k];
				A[ipivot][k]=A[j][k];
				A[j][k]=B[k];
			}
		}
		
		pivot=A[j][j];
		for(k=0;k<nc;k++){
			A[j][k]=A[j][k]/pivot;
		}
		for(i=0;i<nr;i++){
			if(i!=j){
				pivot=A[i][j];
				for(k=0;k<nc;k++){
					A[i][k]-=A[j][k]*pivot;
				}
			}
		}
	}
	
	for(i=0; i<nr; i++)
	{
/*			printf("%d\t %lf\t %lf\n", i, a[i], A[i][nc-1]);
	scanf("%s\n", buffer);*/
	a[i]=A[i][nc-1];
//			printf("%d\t %lf\t %lf\n", i, a[i], A[i][nc-1]);
	}
			
	return 0;
}		


int main(int argc, char *argv[])
{
	FILE *fp;
	FILE *fp1;
	
	char filenamein[100];
	char filenameout[100], filenameout1[100], filenameout2[100], filenameout3[100];
	char buffer[100], buffer1[100], buffer2[100];
	char *path="./LSMO/";
	char *Tstr[19]={"20", "40", "60", "80", "100", "120", "140", "160", "180", "200", "220",
				  "240", "250", "260", "270", "280", "290", "300", "1"};
	int i, isg, j, k, l, m, n, nn, nr, nc, nw, N, NN;
	double c, ch, d, g, p, pi, sh, th, z, z1, z2, z3, zs, zs1, zsmin;
	double ac1, c1, c2, c3, c4, phi, s1, s2, s3, s4, calpha, salpha, s2alpha, c2alpha;
	double Fa[NC], da[NR], b[NR], alpha, beta;
	double x[100], y[100], y1[100], y2[100];
	double T[300], B[300], rhoH[300];
	double temp, field, rho;
	double thickness;
	
	da[0]=0.1;
	da[1]=0.01;
	da[2]=0.001;
	da[3]=0.0001;
	
	pi=M_PI;
	c=pi/180.0;
//	thickness=3300e-8; // [cm] LSMO film thickness
	thickness=3750e-8; // [cm] LCMO film thickness
	
	j=0;
	while((buffer1[j++]=argv[1][j])!='-');
	buffer1[j-1]=0;
	strcpy(filenameout, "least-sq-out-");
	strcat(filenameout, buffer1);
	strcat(filenameout, ".txt");
	fp=fopen(filenameout, "w");
//	printf("%d\t %s\t %s\n", j, argv[1], filenameout);
//	scanf("%s\n", buffer);
	NN=argc-1;	
	for(nn=0;nn<NN;nn++)
	{
		/*  file names from command line  */
		nr=5;
		nc=6;
		j=0;
		while((filenamein[j++]=*(argv[1+nn]++))!=0);
		j=0;
		while((buffer2[j++]=filenamein[j])!='-');
		buffer2[j-1]=0;
		if(strcmp(buffer1, buffer2)!=0)
		{
			fclose(fp);
			strcpy(filenameout, "least-sq-out-");
			strcat(filenameout, buffer2);
			strcat(filenameout, ".txt");
			strcpy(buffer1, buffer2);
			fp=fopen(filenameout, "w");
		}
//		printf("%s\t %s\t %d\n", buffer1, buffer2, strcmp(buffer1, buffer2));
//		scanf("%s\n", buffer);
		
//		printf("%d\t %s\n", nn, filenamein);
//		scanf("%s\n", buffer);
		fprintf(fp, "%s\n", filenamein);
		/*  file names for output */
		strcpy(filenameout1, "fit-");
		strcpy(filenameout2, "para-");
		strcpy(filenameout3, "devi-");
		strcat(filenameout1, filenamein); 
		strcat(filenameout2, filenamein); 
		strcat(filenameout3, filenamein);
		/*  extract values for temperature and field from a filename  */
		i=0;
		while((buffer[i]=filenamein[i++])!='K');
		buffer[i-1]=0;
//		printf("T %s  %ld\n", buffer, strlen(buffer));
//		printf("T %lf\n", atof(buffer));
		T[nn]=atof(buffer);
		j=0;
		while((buffer[j]=filenamein[1+i+j++])!='T');
		buffer[j-1]=0;
		printf("B %s\n", buffer);
		printf("T %lf\n", atof(buffer));
		B[nn]=atof(buffer);
		
		
		/*  read data from an input file  */
		fp1=fopen(filenamein, "r");
		i=0;
		z=0.0;
		while(feof(fp1)==0)
//		while(z<180)
		{
			fscanf(fp1, "%lf\t %lf\t %lf\t %lf\n", &x[i], &y[i], &y1[i], &y2[i]);
			z=x[i];
			if((z>0 && z<179) || (z>181 && z<360))
//			if(z>181)
			{
				printf("%lf\t %lf\t %lf\t %lf\n", x[i], y[i], y1[i], y2[i]);
				i++;
			}
		}
		printf("\n");
		fclose(fp1);
		N=i-1;
		
/*  least square fitting (perpendicular method)  */
/*  set initial approximate values for the parameters a[i]  */
		
		nr=6;
		nc=7;
		a0[0]=0.0015;
		a0[1]=0.0025;
		a0[2]=-0.045;
		a0[3]=0.00004;
		a0[4]=0.035;
		a0[5]=0.55;
		
		zs=0;
		for(i=0;i<N;i++)
		{
			th=x[i]*c;
			if(th>pi)
			{
				sh=-pi;
				ch=-1.0;
			}
			else
			{
				sh=0;
				ch=1.0;
			}
			phi=aniso_magn_r(a0[5], th)+sh;
			c1=cos(th);
			s1=sin(th);
			calpha=cos(th+phi);
			salpha=sin(th+phi);
			c2alpha=cos(2*(th+phi));
			s2alpha=sin(2*(th+phi));
			z=cos(phi)-2*ch*a0[5]*c2alpha;
			z1=s2alpha/z;
			z2=cos(phi)/z;
			alpha=-c*a0[0]*s1-c*(2*a0[1]*s2alpha+a0[2]*salpha)*z2+a0[3];
			beta=-1;
			p=1.0/(alpha*alpha+beta*beta);
			
			Fa[6]=a0[0]*c1+a0[1]*c2alpha+a0[2]*calpha+a0[3]*x[i]+a0[4]-y[i];
			zs+=p*Fa[6]*Fa[6];
		}
		printf("L %lf\n", zs);
//		scanf("%s\n", buffer);
		zs1=zs;
		
		k=0;
		while(k<2)
		{
			for(j=0;j<nr;j++)
			{
				n=0;
				isg=1;
				m=0;
				l=0;
				while(1)
				{	// A
					a0[j]+=isg*da[n];
//					printf("%d\t %d\t %d\t %d\t %lf\n", j, m, n, isg, a0[j]);
					// B
					zs=0;
					sh=0;
					for(i=0;i<N;i++)
					{
						th=x[i]*c;
						if(th>pi)
						{
							sh=-pi;
							ch=-1.0;
						}
						else
						{
							sh=0;
							ch=1.0;
						}
						phi=aniso_magn_r(a0[5], th)+sh;
						c1=cos(th);
						s1=sin(th);
						calpha=cos(th+phi);
						salpha=sin(th+phi);
						c2alpha=cos(2*(th+phi));
						s2alpha=sin(2*(th+phi));
						z=cos(phi)-2*ch*a0[5]*c2alpha;
						z1=s2alpha/z;
						z2=cos(phi)/z;
						alpha=-c*a0[0]*s1-c*(2*a0[1]*s2alpha+a0[2]*salpha)*z2+a0[3];
						beta=-1;
						p=1.0/(alpha*alpha+beta*beta);
						
						Fa[6]=a0[0]*c1+a0[1]*c2alpha+a0[2]*calpha+a0[3]*x[i]+a0[4]-y[i];
						zs+=p*Fa[6]*Fa[6];
					}
/*					printf("%d\t %d\t %d\t %d\t %lf\t L %lf\t L1 %lf\t dL %lf\n", 
						   j, m, n, isg, a0[j], zs, zs1, zs-zs1);
					scanf("%s\n", buffer);*/
					// C
					if(zs<zs1) 
					{	// D
						m=1;
						z1=z;
						l++;
						if(l>10 && n==0) break;
					}
					else
					{	// E
						a0[j]-=isg*da[n];
/*					printf("%d\t %d\t %d\t %d\t %lf Reset\n", j, m, n, isg, a0[j]);
					scanf("%s\n", buffer);*/
						// F
						if(m==0)
						{
							// K
							if(isg<0)
							{
								// L
								isg=1;
								// H
								n++;
								// I
								if(n>3) break;
							}
							// M
							else
							{
								isg=-1;
							}
						}
						// J
						else
						// G
						{	
							// H
							n++;
							// I
							if(n>3) break;
						}
						// J
					}
				}
				
				
			}
			k++;
		}
		printf("L %lf\n", zs);
		zs1=zs;
			printf("a0[i]\n");
			for(i=0;i<nr;i++) printf("%lf\t", a0[i]);
			printf("\n");
		
/* End of the determination of the initial approximate values for a0[i] */		
		
		for(j=0;j<nr;j++) b[j]=a0[j];
		zsmin=zs1;
		m=0;
		l=0;
		nw=0;
		while(nw<201)  /*  repeat 5 times  */
//		while(nw<0) 
		{
			for(i=0;i<nr;i++)
			{
				for(j=0;j<nc;j++)
				{
					A[i][j]=0;
				}
			}
			
			//printf("%d\n", N);
			zs=0;
			sh=0;
			for(i=0;i<N;i++)
			{
				th=x[i]*c;
				sh=-pi*(th>=pi);
				ch=2*(th<pi)-1;
				c1=cos(th);
				s1=sin(th);
				phi=aniso_magn_r(a0[5], th)+sh;
				calpha=cos(th+phi);
				salpha=sin(th+phi);
				c2alpha=cos(2*(th+phi));
				s2alpha=sin(2*(th+phi));
				
				z=cos(phi)-2*ch*a0[5]*c2alpha;
				z1=s2alpha/z;
				z2=cos(phi)/z;
				alpha=-c*a0[0]*s1-c*(2*a0[1]*s2alpha+a0[2]*salpha)*z2+a0[3];
				beta=-1;
				p=1.0/(alpha*alpha+beta*beta);
				Fa[0]=c1;
				Fa[1]=c2alpha;
				Fa[2]=calpha;
				Fa[3]=x[i];
				Fa[4]=1.0;
				Fa[5]=-ch*(2*a0[1]*s2alpha+a0[2]*salpha)*z1;
				Fa[6]=a0[0]*c1+a0[1]*c2alpha+a0[2]*calpha+a0[3]*x[i]+a0[4]-y[i];
				Fa[6]*=-1;
				
				zs+=p*Fa[6]*Fa[6];
				
				for(j=0;j<nr;j++)
				{
					for(k=j;k<nc;k++) A[j][k]+=p*Fa[j]*Fa[k];
				}
			
/*			
			for(j=0;j<nr;j++) printf("%lf\t", Fa[j]); printf("\n");
			printf("%lf\t %lf\t %lf\n", x[i], th, phi);
			for(j=0;j<nr;j++)
			{
				for(k=0;k<nc;k++)
				{
					printf("%lf\t ", A[j][k]);
				}
				printf("\n");
			}
			printf("\n");
			scanf("%s\n", buffer);*/
			}
			
			for(j=1;j<nr;j++)
			{
				for(k=0;k<j;k++)
				{
					A[j][k]=A[k][j];
				}
			}
			
/*			printf("%lf\t %lf\t %lf\n", x[i], th, phi);
			for(i=0;i<nr;i++)
			{
				for(j=0;j<nc;j++)
				{
					printf("%lf\t ", A[i][j]);
				}
				printf("\n");
			}
			printf("\n");
			scanf("%s\n", buffer);*/
			
			n=LeastSquare(nr);

			printf("a0[i]\n");
			for(i=0;i<nr;i++) printf("%lf\t", a[i]);
			printf("\n");
			for(i=0;i<nr;i++) printf("%lf\t", a0[i]);
			printf("\n");
			for(i=0;i<nr;i++) a0[i]+=0.02*a[i];
			
			for(i=0;i<nr;i++) printf("%lf\t", a0[i]);
			printf("\n");
			for(i=0;i<nr;i++) printf("%lf\t", b[i]);
			printf("\n");
			printf("m %d\t l %d\t L %lf\t L1 %lf\t dL %lf\t zsmin %lf\n", m, l, zs, zs1, zs-zs1, zsmin);
/*	*/		zs1=zs;
			if(zs1<zsmin) zsmin=zs1;
				
/*			printf("\n%d debug 3 ", nw);
			scanf("%s\n", buffer);*/
			
/*			if(zs<zs1)
			{
				m++;
				zs1=zs;
				if(zs1 < zsmin)
				{
					zsmin=zs1;
					for(j=0;j<nr;j++) b[j]=a0[j];
				}
			}
			else
			{
				l++;
				if(l>150)
				{
					break;
				}
			}*/
			
/*	*/		if(a0[5]<-0.48) a0[5]=-0.48;
			if(a0[5]>1.5) a0[5]=1.5;
			
			nw++;
		}
/*			for(i=0;i<nr;i++) printf("%lf\t", b[i]);
			printf("\n");
			for(j=0;j<nr;j++) a0[j]=b[j];*/
		
		phi=aniso_magn_r(a0[5], 0);
		rhoH[nn]=a0[0]*thickness*1e6;  				// [micro ohm cm]
//		rhoH[nn]=(a0[0]+a0[2]*cos(phi))*thickness*1e6;  				// [micro ohm cm]
//		printf("%lf\t %lf\t %lf\n", T[nn], B[nn], rhoH[nn]);
		
		for(i=0; i<nr; i++)
		{
			fprintf(fp, "%lf", a0[i]);
			if(i!=nr-1) fprintf(fp, "\t "); 
			else fprintf(fp, "\n");
		}
		for(i=0; i<nr; i++)
		{
			printf("%lf", a0[i]);
			if(i!=nr-1) printf("\t "); 
			else printf("\n");
		}
		fprintf(fp, "%lf\t %lf\t %lf\n", T[nn], B[nn], rhoH[nn]);
		printf("\n%lf\t %lf\t %lf\n", T[nn], B[nn], rhoH[nn]);
		
		/* calculation of dispersion */
		zs=0;
		for(i=0;i<N;i++)
		{
//			th=(x[i]-a0[6])*c;
			th=x[i]*c;
			sh=-pi*(th>pi);
			ch=2*(th<pi)-1;
			phi=aniso_magn_r(a0[5], th)+sh;
			c1=cos(th);
			c2=cos(2*th);
			s1=sin(th);
			calpha=cos(th+phi);
			salpha=sin(th+phi);
			c2alpha=cos(2*(th+phi));
			s2alpha=sin(2*(th+phi));
			z=s2alpha/(cos(phi)-2*ch*a0[5]*c2alpha);
			z1=cos(phi)/(cos(phi)-2*a0[5]*c2alpha);
			alpha=-c*a0[0]*s1-c*(2*a0[1]*s2alpha+a0[2]*salpha)*z2+a0[3];
			beta=-1;
			p=1.0/(alpha*alpha+beta*beta);
			Fa[6]=a0[0]*c1+a0[1]*c2alpha+a0[2]*calpha+a0[3]*x[i]+a0[4]-y[i];
			zs+=p*Fa[6]*Fa[6];
		}
		zs/=N;	/*  dispersion  */
		printf("L %lf\t L1 %lf\t dL %lf\n", zs, zs1, zs-zs1);
			
		fp1=fopen(filenameout2, "w");	/* para- */
		for(i=0; i<nr; i++)
		{
			fprintf(fp1, "%lf", a0[i]);
			if(i!=nr-1) fprintf(fp1, "\t ");
			else fprintf(fp1, "\n");
		}
		fprintf(fp1, "%le\n", zs);
		fclose(fp1);
		
		fp1=fopen(filenameout1, "w");	/* fit- */
		//for(i=0;i<N;i++)
		j=0;
		for(i=0;i<361;i+=2)
		{
			th=i*c;
			if(th>=pi && j==0)
			{
				fprintf(fp1, "\n");
				j=1;
			}
			sh=-pi*(th>=pi);
			ch=2*(th<pi)-1;
			phi=aniso_magn_r(a0[5], th)+sh;
			c1=cos(th);
			calpha=cos(th+phi);
			c2alpha=cos(2*(th+phi));
			z=a0[0]*c1+a0[1]*c2alpha+a0[2]*calpha+a0[3]*i+a0[4];
//			printf("%d\t %lf\t %lf\n", i, phi, z);
			fprintf(fp1, "%d\t %lf\n", i, z);
		}
		fclose(fp1);
		
		fp1=fopen(filenameout3, "w");	/* devi */
		for(i=0;i<N;i++)
		{
			phi=aniso_magn_r(a0[5], th);
			c1=cos(th);
			phi=aniso_magn_r(a0[5], th);
			calpha=cos(th+phi);
			c2alpha=cos(2*(th+phi));
			z=a0[0]*c1+a0[1]*c2alpha+a0[2]*calpha+a0[3]*x[i]+a0[4];
			z1=y[i]-(a0[1]*c2alpha+a0[2]*calpha+a0[3]*x[i]+a0[4]);
			z2=y[i]-(a0[0]*c1+a0[1]*c2alpha+a0[3]*x[i]+a0[4]);
			fprintf(fp1, "%lf\t %lf\t %lf\t %lf\t %lf\t %lf\n", x[i], y[i], z, y[i]-z, z1, z2);
		}
		fclose(fp1);
	}
	
//	printf("debug 4\n");
//	scanf("%s\n", buffer);
	for(i=0;i<argc-1;i++)
	{
		for(j=0;j<i;j++)
		{
			z=T[i];
			z1=B[i];
			z2=rhoH[i];
			if(T[i]<T[j])
			{
				for(k=i;k>j;k--)
				{
					T[k]=T[k-1];
					B[k]=B[k-1];
					rhoH[k]=rhoH[k-1];
				}
				T[j]=z;
				B[j]=z1;
				rhoH[j]=z2;
				j=i;
			}
			else if(T[i]==T[j])
			{
//				printf("%d\t %d\t %lf\n", i, j, T[i]);
				l=0;
				while(T[i]==T[j+l] && j+l < i)
				{
					l++;
				}
//				printf("l %d\n", l);
				for(m=0;m<l;m++)
				{
					if(B[i]<B[j+m] || B[i]==B[j+m])
					{
						for(k=i;k>j+m;k--)
						{
							T[k]=T[k-1];
							B[k]=B[k-1];
							rhoH[k]=rhoH[k-1];
						}
						T[j+m]=z;
						B[j+m]=z1;
						rhoH[j+m]=z2;
						m=l;
					}
					else
					{
						for(k=i;k>j+l;k--)
						{
							T[k]=T[k-1];
							B[k]=B[k-1];
							rhoH[k]=rhoH[k-1];
						}
						T[j+l]=z;
						B[j+l]=z1;
						rhoH[j+l]=z2;
						m=l;
					}
				}
				j=i;
			}
		}
	}
	fclose(fp);
	
/*	for(i=0;i<argc-1;i++)
	{
		printf("%d\t %lf\t %lf\t %lf\n", i, T[i], B[i], rhoH[i]);
	}*/
	i=0;
	while(i<argc-1)
	{
//		printf("\n%s\n", filenamein);
		j=0;
		while(T[i]!=atof(Tstr[j++]));
//		printf("%d\t %s\n", i, Tstr[j-1]);
		strcpy(filenameout, "rhoH-B-");
		strcat(filenameout, Tstr[j-1]);
		strcat(filenameout, "K-Ca.txt");
		fp=fopen(filenameout, "w");
		j=i;
		while(T[i]==T[j])
		{
			printf("%lf\t %lf\t %lf\n", T[j], B[j], rhoH[j]);
			fprintf(fp, "%lf\t %lf\n", B[j], rhoH[j]);
			//scanf("%s\n", buffer);
			j++;
		}
//		printf("close %s\n", filenameout);
		fclose(fp);
		//scanf("%s\n", buffer);
		i=j;
	}
}