clear;

 # --------------------
 # Graphene Properties:
 # --------------------

 graph_name = "C (graphene) - broadband - top";
 graph_side_name = "C (graphene) - broadband - sidewall";
 # scattering rate
 graph_scattering = 0.015;
 # temperature
 graph_temperature = 300;
 # import chemical potential results
 filename_top = "mu_vs_voltage_top.txt";
 V_top = readdata(filename_top);
 filename_side = "mu_vs_voltage_leftwall.txt";
 V_side = readdata(filename_side);
 V_gate_array = pinch(V_top,2,1); 
 npoints = length(V_gate_array);
 myuc_top = pinch(V_top,2,2);
 myuc_side = pinch(V_side,2,2);

 # -------------------
 # Results to Collect:
 # -------------------

 loss_array = matrix(npoints,1);

 # -------------------------
 # Transmission Calculation:
 # -------------------------

 load("graphene_electro-optic_modulator.lms");
 switchtolayout;
 # select FDE solver
 setactivesolver('FDE');
 # set FDE initial guess
 setnamed("FDE","use max index",false);
 setnamed("FDE", "n",1.8);
 setnamed("FDE", "number of trial modes", 1);
 # set constant graphene parameters
 setmaterial(graph_name,"scattering rate (eV)",graph_scattering);
 setmaterial(graph_name,"temperature (K)",graph_temperature);
 setmaterial(graph_side_name,"scattering rate (eV)",graph_scattering);
 setmaterial(graph_side_name,"temperature (K)",graph_temperature);

 for(I = 1:npoints){

    switchtolayout;
     
    # set chemical potential
    setmaterial(graph_name,"chemical potential (eV)",abs(myuc_top(I)));
    setmaterial(graph_side_name,"chemical potential (eV)",abs(myuc_side(I)));
   
    # import charge density profile 
    if ((myuc_top(I)<0.05) & (myuc_side(I)<0.05)){
        ?msg = "using model 1";
        setnamed("np density 1","enabled",true);
        setnamed("np density 1","V_gate",V_gate_array(I));
        setnamed("np density 2","enabled",false);
    }else{
        if ((myuc_top(I)>0.05) & (myuc_side(I)>0.05)){
            ?msg = "using model 2";
            setnamed("np density 1","enabled",false);
            setnamed("np density 2","enabled",true);
            setnamed("np density 2","V_gate",V_gate_array(I));
        }else{
            ?msg = "Warning: two different models might be necessary for the chemical potentials on top and side";
            break;
        }
    }
    setnamed("Graphene-Layer-Top","material",graph_name);
    setnamed("Graphene-Layer-Sidewall","material",graph_side_name);
    
    # run FDE simulation
    findmodes;
    # collect loss
    loss_array(I) = getdata("FDE::data::mode1","loss");
    ?msg = "Loss: " + num2str(loss_array(I));

 } 

 # -----------------------------------
 # Transmission vs. Gate Voltage Plot:
 # -----------------------------------

 # experimental data
 Texp_data = readdata("graphene_electro-optic_modulator_paper_data.txt"); 
 voltage_shift = 0;
 plot(V_gate_array+voltage_shift, -(loss_array-min(loss_array))*1.0e-6, "bias voltage (V)", "transmission (dB/um)");
 holdon;
 plot(pinch(Texp_data,2,1), pinch(Texp_data,2,2),"bias voltage (V)", "transmission (dB/um)"," ","plot points");
 legend("numerical","experimental");
 setplot("x min",-5.5);
 setplot("x max",5);
 holdoff;