//#define debug
#include "thinfilm_BLS.h"

extern int ramp_start;

#ifdef USEDEGREES
#define rtod(x) (x)
#else
#define rtod(x) ((x)/M_PI*180)
#endif

void printmodeprofiles(tf_Layer *l, tf_Global g, double *Freq, int maxidx)
{
    int k,j, wmin;
    double **A,**v,*w,*m, norm_m;
    double z;
    FILE *output;
    double Int[MAXFREQS];

#ifdef mokesim
     if(par.logfile != NULL)
         output = par.logfile;
     else
         output = stdout;
#else
     if(g.outp != NULL)
         output = g.outp;
     else
         output = stdout;
#endif

A = matrix(1,2*g.n_layers,1,2*g.n_layers);
v = matrix(1,2*g.n_layers,1,2*g.n_layers);
w = vector(1,2*g.n_layers);
m = vector(1,2*g.n_layers);

fprintf (output, "print \"\\nB=%.4f\\n\"\n",g.B);

//plot moke-profile
if(g.output_type == 3 && ramp_start == 1)  {

fprintf (output,"set title \"MO r_ps depthprofiles\"\n");
fprintf (output,
"plot \
\"-\" tit \"re_L\" wi lin, \
\"-\" tit \"im_L\" wi lin, \
\"-\" tit \"re_P\" wi lin, \
\"-\" tit \"im_P\" wi lin, \
0 tit \"\" wi lin;");

//real part of long. MOKE
fprintf (output, "\n#\n#\n#y(m)\tre_L(rad/m)\n");
fprintf (output,"%e\t%.5f\n",0.0,0.0);
for (k=1,z=0; k<=g.n_layers; k++)
        {
        fprintf (output,"%e\t%e\n",z,l[k].re_L);
        z += l[k].d;
        fprintf (output,"%e\t%e\n",z,l[k].re_L);
        if(l[k].s != 0.0)
        fprintf (output,"%e\t%e\n",z,0.0);
        z += l[k].s;
        if(l[k].s != 0.0)
        fprintf (output,"%e\t%e\n",z,0.0);
        }
fprintf (output,"%e\t%e\n",z,0.0);
fprintf (output, "e\n#\n#\n");

//imaginary part of long. MOKE
fprintf (output, "\n#\n#\n#y(m)\tim_L(rad/m)\n");
fprintf (output,"%e\t%.5f\n",0.0,0.0);
for (k=1,z=0; k<=g.n_layers; k++)
        {
        fprintf (output,"%e\t%e\n",z,l[k].im_L);
        z += l[k].d;
        fprintf (output,"%e\t%e\n",z,l[k].im_L);
        if(l[k].s != 0.0)
        fprintf (output,"%e\t%e\n",z,0.0);
        z += l[k].s;
        if(l[k].s != 0.0)
        fprintf (output,"%e\t%e\n",z,0.0);
        }
fprintf (output,"%e\t%e\n",z,0.0);
fprintf (output, "e\n#\n#\n");

//real part of polar MOKE
fprintf (output, "\n#\n#\n#y(m)\tre_P(rad/m)\n");
fprintf (output,"%e\t%.5f\n",0.0,0.0);
for (k=1,z=0; k<=g.n_layers; k++)
        {
        fprintf (output,"%e\t%e\n",z,l[k].re_P);
        z += l[k].d;
        fprintf (output,"%e\t%e\n",z,l[k].re_P);
        if(l[k].s != 0.0)
        fprintf (output,"%e\t%e\n",z,0.0);
        z += l[k].s;
        if(l[k].s != 0.0)
        fprintf (output,"%e\t%e\n",z,0.0);
        }
fprintf (output,"%e\t%e\n",z,0.0);
fprintf (output, "e\n#\n#\n");

//imaginary part of polar MOKE
fprintf (output, "\n#\n#\n#y(m)\tim_P(rad/m)\n");
fprintf (output,"%e\t%.5f\n",0.0,0.0);
for (k=1,z=0; k<=g.n_layers; k++)
        {
        fprintf (output,"%e\t%e\n",z,l[k].im_P);
        z += l[k].d;
        fprintf (output,"%e\t%e\n",z,l[k].im_P);
        if(l[k].s != 0.0)
        fprintf (output,"%e\t%e\n",z,0.0);
        z += l[k].s;
        if(l[k].s != 0.0)
        fprintf (output,"%e\t%e\n",z,0.0);
        }
fprintf (output,"%e\t%e\n",z,0.0);
fprintf (output, "e\n#\n#\n");

fprintf (output, "print \"\\npress any to continue\"\npause (-1)\n");
}	//end moke-profiles

if(g.output_type <= 3)	{

fprintf (output,"set title \"in-plane magnetization angle\"\n");
fprintf (output,"plot \"-\" tit \"\" wi lin");
fprintf (output, "\n#\n#\n#y(m)\tin-plane magnetization angle(deg)\n");
fprintf (output,"%e\t%.5f\n",0.0,0.0);
for (k=1,z=0; k<=g.n_layers; k++)
        {
        fprintf (output,"%e\t%.5f\n",z,rtod(l[k].theta_M));
        z += l[k].d;
        fprintf (output,"%e\t%.5f\n",z,rtod(l[k].theta_M));
        if(l[k].s != 0.0)
        fprintf (output,"%e\t%.5f\n",z,0.0);
        z += l[k].s;
        if(l[k].s != 0.0)
        fprintf (output,"%e\t%.5f\n",z,0.0);
        }
fprintf (output,"%e\t%.5f\n",z,0.0);
fprintf (output, "e\n#\n#\n");
fprintf (output, "print \"\\npress any to continue\"\npause (-1)\n");
}	//end outp-type <= 3


for(j=0; j<=maxidx; j++)
{
g.omega = Freq[j];
tf_buildmatrix (A, l, g);
svdcmp (A, 2*g.n_layers, 2*g.n_layers, w, v);
wmin=1;
for (k=1; k<=2*g.n_layers; k++)
        if (fabs (w[k]) < fabs (w[wmin])) wmin=k;

for (k=1; k<=2*g.n_layers; k++)
        m[k] = v[k][wmin];

if(g.output_type == 1) {
norm_m =0;
for (k=1; k<=g.n_layers; k++)
        { norm_m += fabs(m[2*k-1]/l[k].M);
          norm_m += fabs(m[2*k]/l[k].M); }

for (k=1; k<=2*g.n_layers; k++)
        m[k] = m[k]/norm_m;
}

if(g.output_type >= 2)
//normalize modes to energy 1J/m^2
{
for (k=1; k<=g.n_layers; k++) {
        m[2*k-1] = m[2*k-1] /l[k].d / l[k].M;
	m[2*k] = m[2*k] /l[k].d / l[k].M;
}

#ifdef debug
printf("\nFrequency: %g\n",Freq[j]);
for (k=1; k<=2*g.n_layers; k++)
	printf("\n%e ",m[k]);
#endif

norm_m = modeenergy(l, g, A, m);

#ifdef debug
printf("mode#%d-energy %e\n", j, norm_m);
#endif
if(norm_m >= 0)
	norm_m = 1/sqrt(norm_m);
else
	norm_m = 0;

for (k=1; k<=2*g.n_layers; k++)
        m[k] = norm_m * m[k];
}

//now printout solution:
//first mode profiles:
if(g.output_type <= 3)	{

fprintf (output,"set title \"mode#%d profile, %.4fT, %.3fGHz\n",j+1,g.B,Freq[j]);
fprintf (output,
"plot \"-\" tit \"in-plane\" wi lin,\"\" tit \"out-of-plane\" wi lin,0 tit \"\" wi lin\n");

if(g.output_type == 1)
fprintf (output, "#y(m)\tmx/(dMs):\n");
if(g.output_type >= 2)
fprintf (output, "#y(m)\tmx(rad/sqrt(J/m^2)):\n");

fprintf (output,"%e\t%.5f\n",0.0,0.0);
for (k=1,z=0; k<=g.n_layers; k+=1)
{
if(g.output_type == 1)
fprintf (output,"%e\t%.5f\n",z,m[2*k-1]/l[k].d);
if(g.output_type >= 2)
fprintf (output,"%e\t%.5f\n",z,m[2*k-1]);

z += l[k].d;

if(g.output_type == 1)
fprintf (output,"%e\t%.5f\n",z,m[2*k-1]/l[k].d);
if(g.output_type >= 2)
fprintf (output,"%e\t%.5f\n",z,m[2*k-1]);

if(l[k].s != 0.0)
fprintf (output,"%e\t%.5f\n",z,0.0);
z += l[k].s;
if(l[k].s != 0.0)
fprintf (output,"%e\t%.5f\n",z,0.0);
}
fprintf (output,"%e\t%.5f\n",z,0.0);
fprintf (output,"e\n");

if(g.output_type == 1)
fprintf (output, "\n#y(m)\tmy/(dMs):\n");
if(g.output_type >= 2)
fprintf (output, "\n#y(m)\tmy(rad/sqrt(J/m^2)):\n");

fprintf (output,"%e\t%.5f\n",0.0,0.0);
for (k=1,z=0; k<=g.n_layers; k+=1)
{
if(g.output_type == 1)
fprintf (output,"%e\t%.5f\n",z,m[2*k]/l[k].d);
if(g.output_type >= 2)
fprintf (output,"%e\t%.5f\n",z,m[2*k]);

z += l[k].d;

if(g.output_type == 1)
fprintf (output,"%e\t%.5f\n",z,m[2*k]/l[k].d);
if(g.output_type >= 2)
fprintf (output,"%e\t%.5f\n",z,m[2*k]);

if(l[k].s != 0.0)
fprintf (output,"%e\t%.5f\n",z,0.0);
z += l[k].s;
if(l[k].s != 0.0)
fprintf (output,"%e\t%.5f\n",z,0.0);
}
fprintf (output,"%e\t%.5f\n",z,0.0);
fprintf (output, "e\n\n");
fprintf (output, "print \"\\npress any key to continue\"\npause (-1)\n");
}	//end g.output_type <= 3

if(g.output_type >= 3)  {
double re_rps=0, im_rps=0;

for (k=1; k<=g.n_layers; k++)	{
re_rps += (cos(l[k].theta_M)*l[k].re_L * m[2*k-1]
	- l[k].im_P * m[2*k]) *l[k].d;
im_rps += (cos(l[k].theta_M)*l[k].im_L * m[2*k-1] 
	+ l[k].re_P * m[2*k]) *l[k].d;

if(g.output_type >= 5)	{	
re_rps += sqr(cos(l[k].theta_M))*l[k].re_A*m[2*k-1]*l[k].d
		- cos(l[k].theta_M)*l[k].im_C*m[2*k]*l[k].d;		//polar quadratic effect seems to be small
im_rps += sqr(cos(l[k].theta_M))*l[k].im_A*m[2*k-1]*l[k].d
		+ cos(l[k].theta_M)*l[k].re_C*m[2*k]*l[k].d;		//polar quadratic effect seems to be small
}
}
Int[j] = sqr(re_rps) + sqr(im_rps);
}
}


//plot spectra
if(g.output_type >= 3)  {

if(ramp_start == 1)	{
fprintf (output, "PAUSE = -1\n");
fprintf (output, "W=1; set samples 501\n"); }


for(j=0; j<=maxidx; j++)	
fprintf (output, "s%d(x)=(2*I%d*W%d/pi)/(W%d**2+4*(x-F%d)**2); W%d=W\n", j, j, j, j, j,j);

fprintf (output, "spectra(x)=s0(x)");
for(j=1; j<=maxidx; j++)
fprintf (output, "+s%d(x)",j);
fprintf (output, "\n");

for(j=0; j<=maxidx; j++)
fprintf (output, "I%d=%e; F%d=%e\n", j, Int[j], j, Freq[j]);

fprintf (output,"set title \"spectra, %.4fT \"\n",g.B);
fprintf (output,"plot [%e:%e] spectra(x) wi lin;\n", g.f1, g.f2 );
fprintf (output, "print \"\\npress any key to continue\"\npause (PAUSE)\n");
}	


free_matrix(A,1,2*g.n_layers,1,2*g.n_layers);
free_matrix(v,1,2*g.n_layers,1,2*g.n_layers);
free_vector(w,1,2*g.n_layers);
free_vector(m,1,2*g.n_layers);
}

//calculate mode energy (J/m^2)
double modeenergy(tf_Layer *l, tf_Global g, double **A, double *m)
{
int i, j, n = g.n_layers;
double E = 0.0;

g.omega = 0.0;
tf_buildmatrix (A, l, g);

for(i=1;i<=n;i++)
for(j=1;j<=2*n;j++)
{
E += ( 
	  A[2*i-1][j]*m[2*i]*m[j]  
	+ A[2*i][j]*m[2*i-1]*m[j] 
     )
	*l[i].d *l[i].M ;

#ifdef debug
printf("\n%e ",E);
#endif
}

#ifdef debug
printf("\n\n");
printf("\n%e ",0.25*E);
#endif

return(0.25*E);
}

#define YEHMOKEOUTP "./YEHMOKE.dat"

int get_mocoeff(tf_Layer *l,tf_Global *g)
{
int i;
char buffer[LINELENGTH];
FILE *f;
double D[2], re_L[2], im_L[2], re_P[2], im_P[2], 
	re_A[2], im_A[2], re_C[2], im_C[2];
double d_total;

//fprintf(stderr,"\ncalled get_mocoeff\n");

if(g->output_type >= 3)
sprintf(buffer,"yeh-moke-0.15 %s \
\"output_file %s\" \"output_type BLS_depthprofile_linear\"",
	g->yehmoke_file, YEHMOKEOUTP);
if(g->output_type >= 5)
sprintf(buffer,"yeh-moke-0.15 %s \
\"output_file %s\" \"output_type BLS_depthprofile_quadratic\"",
         g->yehmoke_file, YEHMOKEOUTP);

fprintf(stderr,"\ntrying to execute program:\n%s\n\
this may take a while ...\n",buffer);

i=system(buffer);

if(i==-1)	{
 fprintf(stderr,"\nERROR: could not execute %s\n",buffer); 
#ifdef mokesim
 fprintf(stderr,"\npress any key to continue");
 getkey ();	
#else
exit(1);
#endif
return(1); }


if(i!=0)	{
 fprintf(stderr,"\nERROR: \"%s\" returned nonzero value (%d)\n",buffer,i);
#ifdef mokesim
 fprintf(stderr,"\npress any key to continue");
 getkey ();
#else
exit(1);
#endif
return(1); }
	
f=fopen(YEHMOKEOUTP,"r");

if(f == NULL)	 {
 fprintf(stderr,"\nERROR: could not open file %s\n",
	YEHMOKEOUTP);
#ifdef mokesim
 fprintf(stderr,"\npress any key to continue");
 getkey ();
#else
exit(1);
#endif
return(1); }

//starting to read moke-depthprofile:

for(i=0;i<2;i++)
	D[i]=re_P[i]=im_P[i]=re_L[i]=im_L[i]=re_A[i]=im_A[i]=re_C[i]=im_C[i]=0;

while(1)	{
if(lineget(buffer,f) == 0) return(2);
sscanf(buffer,"%le%le%le%le%le%le%le%le%le",
	&D[0],&re_L[0],&im_L[0],&re_P[0],&im_P[0],
	&re_A[0],&im_A[0],&re_C[0],&im_C[0]);
if(re_L[0] != 0.0) break;
}

if(lineget(buffer,f) == 0)return(2);
sscanf(buffer,"%le%le%le%le%le%le%le%le%le",
	&D[1],&re_L[1],&im_L[1],&re_P[1],&im_P[1],
	&re_A[1],&im_A[1],&re_C[1],&im_C[1]);
d_total = 0.5*l[1].d;


for(i=1;i<=g->n_layers;)	{
if(re_L[1] == 0.0)
{
if(lineget(buffer,f) == 0) break;
        sscanf(buffer,"%le%le%le%le%le%le%le%le%le",
	&D[1],&re_L[1],&im_L[1],&re_P[1],&im_P[1],
	&re_A[1],&im_A[1],&re_C[1],&im_C[1]);
        continue;
}

if(fabs(d_total-D[1]) <= fabs(d_total-D[0])) {
	D[0]=D[1];
	re_L[0]=re_L[1];
	im_L[0]=im_L[1];
	re_P[0]=re_P[1];
        im_P[0]=im_P[1];
	re_A[0]=re_A[1];
        im_A[0]=im_A[1];
	re_C[0]=re_C[1];
        im_C[0]=im_C[1];

	if(lineget(buffer,f) == 0) break;
	sscanf(buffer,"%le%le%le%le%le%le%le%le%le",
	&D[1],&re_L[1],&im_L[1],&re_P[1],&im_P[1],
	&re_A[1],&im_A[1],&re_C[1],&im_C[1]);
	continue;
	}
else	{
//printf("d_total=%e D_MO=%e\n",D[0],d_total);
l[i].re_L = re_L[0];
l[i].im_L = im_L[0];
l[i].re_P = re_P[0];
l[i].im_P = im_P[0];
l[i].re_A = re_A[0];
l[i].im_A = im_A[0];
l[i].re_C = re_C[0];
l[i].im_C = im_C[0];

d_total +=  0.5*l[i].d + l[i].s; 
i ++;
if(i<=g->n_layers) d_total += 0.5*l[i].d;
}
}

for(;i<=g->n_layers;i++)
{
l[i].re_L = re_L[0];
l[i].im_L = im_L[0];
l[i].re_P = re_P[0];
l[i].im_P = im_P[0];
l[i].re_A = re_A[0];
l[i].im_A = im_A[0];
l[i].re_C = re_C[0];
l[i].im_C = im_C[0];

}

return(0);
}

