//VerilogA for laser,eml,veriloga

`include "constants.vams"
`include "disciplines.vams"

`define MAX_ARRAY 1000

module eml(Vin,Vlam,Vsig1,Vsig2,Opow,Ophase,Olam,Gnd);
    parameter real L = 300u from [0:inf);
    parameter real W = 1.5u from [0:inf);
    parameter real D = 0.2u from [0:inf);
    parameter real GAM = 0.3 from [0:1);
    parameter real T = 10n from [0:inf);
    parameter real BEATsp = 5e-5 from [0:inf);
    parameter integer ID = 0 from [0:1];
    parameter real G0 = 25000 from [0:inf);
    parameter real Ntr = 1e24 from [0:inf);
    parameter real EPS = 6e-23 from [0:inf);
    parameter real B = 1e-16 from [0:inf);
    parameter real A = 7.5E-41 from [0:inf);
    parameter real Rl = 1e-8 from [0:1);
    parameter real ALFA0 = 5000 from [0:inf);
    parameter real Q1 = 1.9998 from [0:inf);
    parameter real Q2 = 1.4391E-4 from [0:inf);
    parameter real Rat = 1 from [0:inf);
    parameter real Ng = 3.6 from [0:inf);
    parameter real N0 = 6.5e13 from [0:inf);
    parameter real Vbi = 1.13 from [0:inf);
    parameter real EIT = 2 from [0:inf);
    parameter real Rs = 1E-5 from [0:inf);
    parameter real Cp = 0p from [0:inf);
    parameter real Rd = 1E15 from [0:inf);
    parameter real Csc0 = 0p from [0:inf);
	parameter real Gam_c = 0 from (-inf:inf);
	parameter integer tau_L = 6 from [1:inf);
	parameter integer num_carriler = 2 from [1:inf);
	parameter integer num_voltage = 6 from [1:inf);
	//TODO: file for multinomial coefficients
	parameter string tau_file = "./tau_mat.in";
	parameter string alpha_file = "./alpha_mat.in";
	parameter string gamma_file = "./gamma_mat.in";
	//Use LF!!!

    input Vin,Vlam,Vsig1,Vsig2;
    output Opow,Ophase,Olam;
    inout Gnd;
    electrical Vin,Vlam,Vsig1,Vsig2,Opow,Ophase,Olam,Gnd;
    electrical Ns_eam,Nph_eam;
	electrical Nin1,Ns,Nph;

	real Vcalc;
	real Gamma=0,Alpha=0;
	real tau=0;
	real opow,ophase;
	integer file_tau,file_alpha,file_gamma;
	integer i,j;
	real f=1e-18;
	real readfile=0;
	real tau_array[0:`MAX_ARRAY];
	real alpha_array[0:`MAX_ARRAY];
	real gamma_array[0:`MAX_ARRAY];

    real VT=`P_K*$temperature/`P_Q;
    real Vact;
    real Cg;
    real CPL,Tph,QV,Cph,Rph;
    real Cd,Csc;
    real N,G;
	real opow_laser,ophase_laser;
	
	real UL,UW,UD;
	real UN0,UG0,UNtr,UEPS;
	real UA,UB,UALFA0,ULAMBDA,ULSPEED,UQ2;

    branch (Gnd,Ns_eam) Cs,Is;
	branch (Ns_eam,Gnd) It;
	branch (Gnd,Nph_eam) Cph_eam,Iph;

	branch (Nin1,Gnd) R2,C1,C2,C3,B1,B2,B3,B4;
    branch (Gnd,Ns) B5,B6,R3,C4;
	branch (Gnd,Nph) B7,C5;

	analog initial begin
		UL=L*1e6;
		UW=W*1e6;
		UD=D*1e6;
		UG0=G0*1e-6;
		UN0=N0*1e-18;
		UNtr=Ntr*1e-18;
		UEPS=EPS*1e18;
		UA=A*1e36;
		UB=B*1e18;
		UALFA0=ALFA0*1e-6;
		ULSPEED=`P_C*1e6;
		UQ2=Q2*1e6;

        Vact=UL*UW*UD;
        Cg=ULSPEED/Ng;
        Tph=Ng/(ULSPEED*(UALFA0+Q1/UQ2));
        QV=`P_Q*Vact;
        Cph=QV;
        Rph=Tph/QV;
	end

	analog begin
		ULAMBDA=V(Vlam)*1e6;
		CPL=`P_H*Cg*Vact*Rat*ULSPEED*(1-Rl)/ULAMBDA/UQ2;
 		N=UN0*(exp(V(Nin1,Gnd)/EIT/VT)-1.0);

        if (N<UNtr)
            G=0;
        else begin
            if (ID<1)
                G=UG0*(N/UNtr-1.0)/(1.0+UEPS*abs(V(Ns,Gnd)));
            else
                G=UG0*ln(N/UNtr)/(1.0+UEPS*abs(V(Ns,Gnd)));
        end

        Cd=QV*UN0*exp(V(Nin1,Gnd)/EIT/VT)/(EIT*VT);
        if (V(Nin1,Gnd)<Vbi)
            Csc=Csc0/sqrt(1.0-V(Nin1,Gnd)/Vbi);
        else
            Csc=Csc0/sqrt(0.1);

        V(Vin,Nin1) <+ I(Vin,Nin1)*Rs;
        V(R2) <+ I(R2)*Rd;
        I(C1) <+ Cp*ddt(V(C1));
        I(C2) <+ Cd*ddt(V(C2));
        I(C3) <+ Csc*ddt(V(C3));
        I(B1) <+ QV*N/T;
        I(B2) <+ QV*UA*pow(N,3);
        I(B3) <+ QV*UB*pow(N,2);
        I(B4) <+ QV*Cg*GAM*G*abs(V(Gnd,Ns));

        I(B5) <+ BEATsp*QV*UB*pow(N,2);
        I(B6) <+ QV*Cg*GAM*G*abs(V(Gnd,Ns));
        I(C4) <+ Cph*ddt(V(C4));
        V(R3) <+ I(R3)*Rph;

		//I(B7) <+ te_ph1*ALPHA/2;
		//I(C5) <+ ddt(V(C5));

		opow_laser = CPL*V(Ns,Gnd);
		ophase_laser=0;

		if (readfile==0) begin
			file_tau=$fopen(tau_file,"r");
			file_alpha=$fopen(alpha_file,"r");
			file_gamma=$fopen(gamma_file,"r");
			if (file_tau==0 || file_alpha==0 || file_gamma==0) begin
				$display("Error: Could not open file");
				$finish;
			end
			for (i=0;i<tau_L;i=i+1) begin
				$fscanf(file_tau,"%f",tau_array[i]);
			end
			for (i=0;i<num_voltage;i=i+1) begin
				for (j=0;j<num_carriler;j=j+1) begin
					$fscanf(file_alpha,"%f",alpha_array[i*num_carriler+j]);
					$fscanf(file_gamma,"%f",gamma_array[i*num_carriler+j]);
				end
			end
			$fclose(file_tau);
			$fclose(file_alpha);
			$fclose(file_gamma);
			readfile=1;
		end

		Vcalc=V(Vsig1)-V(Vsig2);
		tau=0;
		for (i=0;i<tau_L;i=i+1) begin
			tau=tau+tau_array[i]*pow(Vcalc,i);
		end
		Gamma=0;
		Alpha=0;
		for (i=0;i<num_voltage;i=i+1) begin
			for (j=0;j<num_carriler;j=j+1) begin
				Gamma=Gamma+gamma_array[i*num_carriler+j]*pow(Vcalc,j)*pow(V(Ns_eam),i); //pow(N,i)*pow(f,i)
				Alpha=Alpha+alpha_array[i*num_carriler+j]*pow(Vcalc,j)*pow(V(Ns_eam),i); //Caution the unit
			end
		end
		if (Gamma==0)
			opow=0;
		else begin
			opow=pow(sqrt(opow_laser)*exp(-Gamma/2),2);
		end
		
		I(Is) <+ f*V(Vlam)/(`P_H*`P_C)*exp(-Gam_c)*(1-exp(-Gamma+2*Gam_c))*opow_laser;
		I(It) <+ V(Ns_eam)/tau;
		I(Cs) <+ ddt(V(Cs));

		I(Iph) <+ Alpha/2*ddt(Gamma);
		I(Cph_eam) <+ ddt(V(Cph_eam));
		ophase=V(Nph_eam)+ophase_laser;

		V(Opow) <+ opow;
		V(Ophase) <+ ophase/(2*`M_PI)*360.0;
		V(Olam) <+ V(Vlam);
	end

endmodule