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

#define with_numrec /*use nr zbrent function for root finding */
//#define without_numrec /* dont use it -- slower */

#ifdef with_numrec

#define Tol_Nu EPS

tf_Layer *glob_l;
tf_Global *glob_g;
double **glob_m;
double det_power;

extern double nr_det(double nu);
#endif  /* with_numrec */


void cp_huge(Huge_real source, Huge_real *dest)
{dest->log=source.log, dest->sign=source.sign; }


int tf_findfreq (tf_Layer *l, tf_Global *g, double *freq)
{
  double **m;
  int i=0, n_mat; 

  n_mat = 2*g->n_layers;
  //printf("\nn_mat=%d\n",n_mat);
  m = matrix(1,n_mat,1,n_mat);


if(g->n_layers==1)
{
double p,q;

i=0;
g->omega = 0;
tf_buildmatrix (m, l, *g);
p = -m[1][1] -m[2][2];
q  = m[1][1] *m[2][2] -m[1][2] *m[2][1];
//printf("\np=%e,q=%e\n",p,q);

freq[i] = l[1].gamma *(-p/2 +sqrt(sqr(p)/4 -q));
if((freq[i] >= g->f1)&&(freq[i] <= g->f2)) i++;

freq[i] = l[1].gamma *( -p/2 -sqrt(sqr(p)/4 -q));
if((freq[i] >= g->f1)&&(freq[i] <= g->f2)) i++;

//printf("\nnu0=%e,nu1=%e\n",freq[0],freq[1]);
//i=2;
}

if(g->n_layers==2)
{
double p[5], omega[5], omegai[5];
int j;
i=0;
g->omega = 0;
tf_buildmatrix (m, l, *g);

p[2] = -m[3][4]*m[4][3]*sqr(l[2].gamma) 
-(m[4][1]*m[1][4] +m[3][1]*m[1][3] +m[3][2]*m[2][3] +m[4][2]*m[2][4])*l[2].gamma*l[1].gamma
-m[2][1]*m[1][2]*sqr(l[1].gamma);

p[1] = sqr(l[1].gamma)*l[2].gamma*(-m[2][1]*m[3][2]*m[1][3] 
	-m[3][1]*m[1][2]*m[2][3] -m[4][1]*m[1][2]*m[2][4] -m[2][1]*m[4][2]*m[1][4])
+ sqr(l[2].gamma)*l[1].gamma*(-m[4][2]*m[2][3]*m[3][4] 
	-m[3][1]*m[1][4]*m[4][3] -m[3][2]*m[2][4]*m[4][3] -m[4][1]*m[1][3]*m[3][4]);
p[0] = sqr(l[1].gamma)*sqr(l[2].gamma)*( -m[3][1]*m[1][2]*m[2][4]*m[4][3]
	-m[4][1]*m[1][2]*m[2][3]*m[3][4] -m[2][1]*m[4][2]*m[1][3]*m[3][4]
	+m[3][1]*m[4][2]*m[1][3]*m[2][4] -m[3][1]*m[4][2]*m[1][4]*m[2][3]
	-m[2][1]*m[3][2]*m[1][4]*m[4][3] -m[4][1]*m[3][2]*m[1][3]*m[2][4]
	+m[2][1]*m[1][2]*m[3][4]*m[4][3] +m[4][1]*m[3][2]*m[1][4]*m[2][3]);

p[3] = 0;
p[4] = 1;

zrhqr (p, 4, omega, omegai);

for (j=1 ;j < 5; j++)
            { freq[i] = omega[j];
	      if((freq[i] >= g->f1)&&(freq[i] <= g->f2)) i++;} 	
}

if(g->n_layers>2)
{
//scan frequency form f1 to f2 in steps
//to find the solutions det=0
  Huge_real det1, det2, det3,olddet2;
  double detmin;
  double f1, f2, f3;
  double oldf2;

#ifdef debug
printf("\n in thinfilm_findfreq\n");
#endif

#define DF 0.00001 	//frequency precision

  i = 0;
  f1 = g->omega = g->f1 + 1E-10;
  tf_buildmatrix (m, l, *g);
  tf_det (m,n_mat,&det1);
#ifdef debug
printf("%e\t%e\t%e\n",g->omega,det1.log,det1.sign);  
#endif
    for (g->omega += g->d_f; g->omega < g->f2; g->omega += g->d_f)
    {
    f2 = g->omega ;
    tf_updateM (m,l,*g);
    tf_det (m,n_mat,&det2);
#ifdef debug
printf("%e\t%e\t%e\n",g->omega,det2.log,det2.sign);
#endif
    if(det1.sign*det2.sign == -1) //sign change of det


#ifdef without_numrec

  {
#ifdef debug
        printf("\t# sign change  %f",f1);
#endif
        //remenber old values
        cp_huge(det2,&olddet2);
        oldf2 = f2;

NEXTLOOP:
        g->omega = f3 = 0.5*(f2+f1);
//	printf("%e\t%e\t%e\n",g->omega,det2.log,det2.sign);
        tf_updateM (m, l, *g);
        tf_det (m,n_mat,&det3);
#ifdef debug
        printf("\n\t\t\t loop:f1 %f f2 %f f3 %f\n",f1,f2,f3);
        printf("\n\t\t\t    sign1 %d sign2 %d sign 3 %d\n",
			det1.sign,det2.sign,det3.sign);   	
#endif
  if((fabs(f1-f2) < DF)||(det3.sign == 0))
           //found minima
	{
            freq[i]=0.5*(f1+f2);
            i++;
            g->omega = f1 = oldf2;
            cp_huge(olddet2,&det1);
            continue;
        }
  else        
        {
        if (det1.sign*det3.sign == -1) //sign change between f1 and f3
                {f2 = f3; cp_huge(det3,&det2); goto NEXTLOOP;}
        else {
              //sign change between f2 and f3
             f1 = f3; cp_huge(det3,&det1); goto NEXTLOOP;
             }
        }
   }
#endif  /* without_numrec */

#ifdef with_numrec
{
#ifdef debug
        printf("\t# sign change  %f",f1);
#endif

glob_l = l;
glob_g = g;
glob_m =m;
det_power = fmax(det1.log, det2.log);

freq[i] = zbrent(nr_det, f1, f2, Tol_Nu);
#ifdef debug
        printf("solution: %f GHz",freq[i]);
#endif
i++;
}
#endif  /* with_numrec */

 f1=f2;  
 cp_huge(det2,&det1);
}//end frequency scan
}//end g->layers != 1

  free_matrix(m,1,n_mat,1,n_mat);
  return (i - 1);
}

#ifdef with_numrec
double nr_det(double nu)
{
Huge_real det;

glob_g->omega = nu;
tf_updateM (glob_m, glob_l, *glob_g);
tf_det (glob_m, 2*glob_g->n_layers, &det);

return (det.sign * exp(det.log - det_power +10));
};

#endif /* with_numrec */
