#include <math.h>
#include "bmsinc.h"

int main()
{ 
double R=5.0;         // Radius of bent waveguide

int layer = 1;        // no. of internal layer
double d = 0.5;       // width of bend waveguide core

double lambda=1.1;    // vaccume wavelength

Dvector Bhx(layer+1); // locations of dielectric layer interfaces
Dvector Bn(layer+2);  //  refractive indices for dielectric layers

// dielectric layer interfaces are defined with respect to radius R

Bhx(0) =  -d;         // inner layer position= R-d
Bhx(1) =  0.0;        // outer layer position= R

Bn(0) = 1.0;         // = n_s (cladding refractive index)
Bn(1) = 1.5;         // = n_f (core refractive index)
Bn(2) = 1.0;         // = n_c (cover refractive index)

// Define bent waveguide
BSlabWaveguide bwg(layer, Bhx, lambda, Bn, R);


// Required parameters for mode solver
BSLAMS_Parameters bpar;

// Define container to hold modes of bent waveguide
BSlabModeArray bma;
        
// finding modes of BWG
    
// Roots of dispersion equation are searched in complex plane
// in region [b0, b1]x[a0, a1].

double b0 = 5.2;  // lower limit for beta
double b1 = 6.5;  // upper limit for beta
    
double a0 = -6e-3; // lower limit for alpha
double a1 = -1e-3; // upper limit for alpha

int Nb = 10;       // [b0, b1] is divided in Nb points.
int Na = 6;        // [a0, a1] is divided in Na points.

// finding modes of bent waveguide in [b0, b1]x[a0, a1]
int Nbm=bend_mode_analysis(bwg, TE, b0, b1, Nb, a0, a1,Na, bpar, bma);
    

// to power normalise modes
double tmp_totpower;
       
for (int j= 0; j < Nbm; j++)
{
  // compute total modal power using analytic method
  // at an angular posotion theta=0.0
  tmp_totpower =bma(j).totpower(0.0); 
 
  bma(j).normalize(tmp_totpower); // to normalize modal power
}


// Display propagation constants
 cout <<"\nPropagation constants of mode(s):\n";
for (int j= 0; j < Nbm; j++)
{
    cout <<"Mode " << j
	 <<"\t gamma=(" <<bma(j).gamma.re <<", " << bma(j).gamma.im <<")"
	 <<"\t beta=" << bma(j).beta
	 <<"\t alpha=" << bma(j).alpha <<"\n";
} 

cout <<"\n";

// Field values of TE0 mode 
double r, p, x, z;
r=5.385164807;      // radial distance
p=0.380506377;   // angular position = pi/2
x=5.0;
z=2.0;

Complex f;

f=bma(0).field_cart(EY, x, z);
cout <<"Complex valued y component of electric field at"
<<" (x, z)= (" << x <<", " << z <<") : (" << f.re 
<<", " << f.im <<") \n";

f=bma(0).field(EY, r, p);
cout <<"Complex valued y component of electric field at"
<<" (r, p)= (" << r <<", " << p <<") : (" << f.re 
<<", " << f.im <<") \n";

cout <<"Real part of y component of electric field at"
<< " r=" << r <<": " << bma(0).field(EY, REP, r) <<"\n";
  
// Interval on which field is plotted
Interval Bdisplay(bwg.R - 4.0, bwg.R+ 5.0);  

for (int j= 0; j < Nbm; j++)
{
  //plot of real, imaginary and absolute value and phase of field
  bma(j).mfile(EY, 1.0, Bdisplay, 600, '0', '0'+j, 'r'); //real+im
  
  //plot of absolute value
  bma(j).mfile(EY, MOD, 1.0, Bdisplay, 1000, '0', '0'+j, 'L'); 
  
}

Cvector amp(1); 
amp(0)=CC1;   // Amplitude of Fundamental mode

double x0, x1, z0, z1; // defining window
int Nx, Nz;  // defining discrtization points

x0= 0.0;
x1= 7.0;
Nx=200;
z0= -7.0;
z1=7.0;
Nz=400;

// plot of absolute value of EY over [x0, x1]x[z0, z1] window
bma.plot(amp, EY, MOD, x0, x1, Nx, z0, z1, Nz, '0', '0');

// plot of real value of EY over [x0, x1]x[z0, z1] window
bma.plot(amp, EY, REP, x0, x1, Nx, z0, z1, Nz, '0', '0');

// animation of propagation of real value of EY over 
// [x0, x1]x[z0, z1] window
bma.movie(amp, x0, x1, Nx, z0, z1, Nz, 40, '0', '0');
}