//#define debug
//this defines interface functions to use the code from
//mokesim
#include "thinfilm_BLS.h"

#ifdef LOCAL
//import homedirectory-name from main 
//if localinstallation
extern char *home;
#endif

//imported from twisted.c
#define TW_NMAX 200

extern int n1,n2,n3,N;    //point index layer1: 0...n1-1
                //              layer2: n1...n2-1
                //N=n1+n2
extern double dd1,dd2,dd3;
extern double tw_oldang[TW_NMAX+1];

extern double frac(double d)
	{ return(d - (int)(d));}

void
mks2tfBLS (Parameter par,tf_Layer *l,tf_Global *g,double angle1,double angle2,double angle3) 
{
// !!!!!!!!!!!!!!!!!!!!!!!!!!!
//rezende uses positive j2 for biquadratic couplig !!
// !!!!!!!!!!!!!!!!!!!!!!!!!!!

int i;

if(par.signalopt == -2) perp_fmr=1;
else perp_fmr=0;

if(par.layers == 1)
g->n_layers = par.vir_n1;
if(par.layers == 2)
g->n_layers = par.vir_n1 + par.vir_n2;
if(par.layers == 3)
g->n_layers = par.vir_n1 + par.vir_n2 + par.vir_n3;
if(g->n_layers>tf_MAXL)
        {printf("error:r\ntoo many layers\n"),
        exit(1);}


g->theta_B = par.fie_ang;
g->B = par.B;
g->q = par.q;
g->theta_q = par.q_ang;
g->omega = 0;
g->f1 = par.frmin;
g->f2 = par.frmax;
g->d_f = (par.frmax-par.frmin)/par.frsteps;


for(i=1;i<= par.vir_n1;i++)
	{
        l[i].s = l[i].j2 = 0;
        l[i].d = par.d1/(par.vir_n1);
        l[i].K = par.K1;
	l[i].U = par.U1;
        l[i].X = par.X1;
        l[i].M = par.M1;
        l[i].D = par.D1;
        l[i].gamma = par.gamma1; 
        l[i].theta_K =  par.ea_Kc1 + par.sa_ang;
	l[i].theta_U =  par.ea_Ku1 + par.sa_ang;
        l[i].theta_X =  par.ebd1 + par.sa_ang; 
//only uniform magnetization at the moment!
        l[i].theta_M = angle1;
        l[i].j1 = par.D1 *par.M1 /l[i].d; 
        l[i].j2 = 0;
	l[i].j3 = 0;
        l[i].S = 0;
	}

if(par.layers > 1)
for(i=par.vir_n1+1;i<= par.vir_n1+par.vir_n2;i++)
        {
        l[i].s = l[i].j2 = 0;
        l[i].d = par.d2/(par.vir_n2);
        l[i].K = par.K2;
	l[i].U = par.U2;
        l[i].X = par.X2;
        l[i].M = par.M2;
        l[i].D = par.D2;
        l[i].gamma = par.gamma2;
        l[i].theta_K = par.ea_Kc2 + par.sa_ang;
	l[i].theta_U = par.ea_Ku2 + par.sa_ang;
        l[i].theta_X =  par.ebd2 + par.sa_ang;
        l[i].theta_M = angle2;
        l[i].j1 = par.D2 * par.M2 /l[i].d;
        l[i].j2 = 0;
        l[i].j3 = 0;
	l[i].S = 0;
        }

if(par.layers > 2)
for(i=par.vir_n1+par.vir_n2+1;i<= g->n_layers;i++)
        {
        l[i].s = l[i].j2 = 0;
        l[i].d = par.d3/(par.vir_n3);
        l[i].K = par.K3;
        l[i].U = par.U3;
        l[i].X = par.X3;
        l[i].M = par.M3;
        l[i].D = par.D3;
        l[i].gamma = par.gamma3;
        l[i].theta_K = par.ea_Kc3 + par.sa_ang;
        l[i].theta_U = par.ea_Ku3 + par.sa_ang;
        l[i].theta_X =  par.ebd3 + par.sa_ang;
        l[i].theta_M = angle3;
        l[i].j1 = par.D3 * par.M3 /l[i].d;
        l[i].j2 = 0;
        l[i].j3 = 0;
        l[i].S = 0;
        }

//now the coupling interface anisotropy
//and spacer thickness per layer

l[1].S += par.Sb1;
l[par.vir_n1].S += par.St1;

if(par.layers >1)
{
l[par.vir_n1].j1 = par.j11;
l[par.vir_n1].j2 = par.j21;
l[par.vir_n1].j3 = 0;

l[par.vir_n1].s  = par.ds1;
l[par.vir_n1 +1].S += par.Sb2;
l[par.vir_n1+par.vir_n2].S += par.St2;
}

if(par.layers > 2)
{
l[par.vir_n1+par.vir_n2].j1 = par.j12;
l[par.vir_n1+par.vir_n2].j2 = par.j22;
l[par.vir_n1+par.vir_n2].j3 = 0;

l[par.vir_n1+par.vir_n2].s  = par.ds2;
l[par.vir_n1+par.vir_n2+1].S += par.Sb3;
l[par.vir_n1+par.vir_n2+par.vir_n3].S += par.St3;
}


if(par.tw_activate >0)
//set up individual angles if
//twisted magnetization is enabled
{
double d;
int j;

#ifdef debug
for(i=1;i<=N;i++)
  fprintf(stderr,"\n%i %e",i,tw_oldang[i]/M_PI*180);
fflush(stderr);
#endif


//layer1
d=0.5*l[1].d;
for( i=1;i<=par.vir_n1;i++)
  {
  //closest sublayer in twisted calc
  j= (int)(d/dd1 +1); 

  l[i].theta_M = //twisting angle in the middle of the layer l[i]
          tw_oldang[j] +frac(d/dd1 +1)*(tw_oldang[j+1] - tw_oldang[j]);    

  d += l[i].d;
  }

//layer2
d=0.5*l[1+par.vir_n1].d;
for( i=1;i<=par.vir_n2;i++)
  {
  //closest sublayer in twisted calc
  j= (int)(d/dd2 +1);

  l[i+par.vir_n1].theta_M = //twisting angle in the middle of the layer l[i]
         tw_oldang[j+n1] +frac(d/dd2 +1)*(tw_oldang[j+n1+1]-tw_oldang[j+n1]); 

  d += l[i+par.vir_n1].d;   
  }

//layer3
d=0.5*l[1+par.vir_n1+par.vir_n2].d;
for( i=1;i<=par.vir_n3;i++)
  {
  //closest sublayer in twisted calc
  j= (int)(d/dd3 +1);

  l[i+par.vir_n1+par.vir_n2].theta_M = //twisting angle in the middle of the layer l[i]
         tw_oldang[j+n1+n2] +frac(d/dd3 +1)*(tw_oldang[j+n1+n2+1]-tw_oldang[j+n1+n2]);

  d += l[i+par.vir_n1+par.vir_n2].d;
  }
}       //twisted mode


#ifdef debug
for(i=1;i<=g->n_layers;i++)
  fprintf(stderr,"\n%i %e",i,l[i].theta_M /M_PI *180);
fflush(stderr);
#endif

}

int thinfilm_BLS	//return value is the number of solutions
 (Parameter par, double *Freq, double angle1, double angle2, double angle3)
{
int i,j;
tf_Layer l[tf_MAXL+1];
tf_Global g;

mks2tfBLS (par,l,&g,angle1,angle2,angle3);
#ifdef debug
showvalues(l, &g);
#endif

i = tf_findfreq (l, &g, Freq);

#ifdef debug
for(j=0;j<=i;j++)
	fprintf(stderr,"\nnu[%d] = %f",j,Freq[j]);
#endif

if(par.blsamp == 1)	//solve for precesion amplitudes
{
    int k,wmin;
    double **A,**v,*w,*m, norm_m;
    double z;
    FILE *output;
     if(par.logfile != NULL)
         output = par.logfile;
     else
         output = stdout;


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);

tf_buildmatrix (A, l, g);

fprintf (output, "print \"\\nB=%.4f\\n\"\n",g.B);
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,l[k].theta_M/M_PI*180);
	z += l[k].d;
	fprintf (output,"%e\t%.5f\n",z,l[k].theta_M/M_PI*180);
	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");


for(j=0; j<=i; 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;

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

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;

//now printout solution:

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");
fprintf (output, "#y(m)\tmx/(dMs):\n");
fprintf (output,"%e\t%.5f\n",0.0,0.0);
for (k=1,z=0; k<=g.n_layers; k+=1)
{
fprintf (output,"%e\t%.5f\n",z,m[2*k-1]/l[k].d);
z += l[k].d;
fprintf (output,"%e\t%.5f\n",z,m[2*k-1]/l[k].d);
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");

fprintf (output, "\n#y(m)\tmy/(dMs):\n");
fprintf (output,"%e\t%.5f\n",0.0,0.0);
for (k=1,z=0; k<=g.n_layers; k+=1)
{
fprintf (output,"%e\t%.5f\n",z,m[2*k]/l[k].d);
z += l[k].d;
fprintf (output,"%e\t%.5f\n",z,m[2*k]/l[k].d);
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");

}

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);
}

fprintf (stderr, "*");
fflush (stderr);
return (i);
}

void  showvalues(tf_Layer *l, tf_Global *g)
{
int i;
        for(i=1;i<=g->n_layers;i++)
        {printf("\nlayer %d\t s=%e d=%e K=%e M=%e\nj1=%e j2=%e gamma=%e"
                ,i,l[i].s,l[i].d,l[i].K, l[i].M,l[i].j1,l[i].j2,l[i].gamma);
         printf(" theta_Kc=%e theta_M=%e D=%e",l[i].theta_K,l[i].theta_M,l[i].D);
        }
printf("i\n\n q=%e B=%e theta_B=%e F=%e\n",g->q,g->B,g->theta_B,g->omega);
}
