//VerilogA for passive,mmi2x2,veriloga
//Xuanqi Chen, Zhifei Wang, Yi-Shing Chang, Jiang Xu, Jun Feng, Peng Yang, Zhehui Wang, Luan H. K. Duong, ”Modeling and Analysis of Optical Modulators Based on Free-Carrier Plasma Dispersion Effect,” IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems (TCAD), vol. 39, no. 5, pp. 977-990, May 2020.

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

`define MAX_MODES 200
`define MAX_ITER 10000
`define IDX_NUM 50000
`define TOLERANCE 0.0000000001

module mmi2x2(Ipow1,Iphase1,Ilam1,Ipow2,Iphase2,Ilam2,Opow1,Ophase1,Olam1,Opow2,Ophase2,Olam2);
    parameter real design_lam = 1.55e-6 from (0:inf);
    parameter real n_ps = 3.45 from (0:inf);
    parameter real n_sub = 1.45 from (0:inf);
    parameter real width_MMI = 3e-6 from (0:inf);
    parameter real width_port = 0.45e-6 from (0:inf);
    parameter real port_offset = 0.75e-6 from (0:inf);
    parameter real thickness = 0.22e-6 from (0:inf);
    parameter real alpha_mmi = 80 from (0:inf);

    input Ipow1,Iphase1,Ilam1,Ipow2,Iphase2,Ilam2;
    output Opow1,Ophase1,Olam1,Opow2,Ophase2,Olam2;
    electrical Ipow1,Iphase1,Ilam1,Ipow2,Iphase2,Ilam2,Opow1,Ophase1,Olam1,Opow2,Ophase2,Olam2;    

    real ulam, ulam0, uwidth_in, uwidth_MMI, uoffset, uthickness, alpha0;
    real k1_ref, kc_ref;
    integer mode_num, mode_num_ref;
    real beta_x_MMI, beta_in, beta_x_in;
    real neff_MMI[0:`MAX_MODES];
    real beta_MMI[0:`MAX_MODES];
    real length_MMI;
    real ky_in, gamma_in;
    real ky_MMI[0:`MAX_MODES];
    real gamma_MMI[0:`MAX_MODES];
    real c0;
    real c0_dis[0:`MAX_MODES];
    real icoff[0:`MAX_MODES];
    real delta_beta;
    real ocoff_r[0:`MAX_MODES];
    real ocoff_i[0:`MAX_MODES];
    real c1r, c1i, c2r, c2i;
    real in1r, in1i, in2r, in2i, out1r, out1i, out2r, out2i;
    real sum_outfield1, sum_outfield2, sum_outfield3, sum_outfield4;
    real sum_singlemode1,sum_singlemode2,sum_singlemode3,sum_singlemode4;

    real border;

    integer i, j;
    real p=1;
    real h;

    analog function real findzeropoint_rect_x;
        input x_lower_input, x_upper_input, Vx;
        real x_lower_input, x_upper_input, Vx;
        real iter_count;
        real x_lower, x_upper;
        real y_mid;

        real y_lower, x_mid;

        begin
            findzeropoint_rect_x = -1;
            x_lower = x_lower_input;
            x_upper = x_upper_input;
            iter_count = 0;
            y_mid = 1;
            while (iter_count<`MAX_ITER && abs(y_mid) > `TOLERANCE) begin
                x_mid = (x_lower + x_upper) / 2;
                y_lower = tan(x_lower / 2) - sqrt(pow(Vx,2) - pow(x_lower,2)) / x_lower;
                y_mid = tan(x_mid / 2) - sqrt(pow(Vx,2) - pow(x_mid,2)) / x_mid;

                if (abs(y_mid) < `TOLERANCE) begin
                    findzeropoint_rect_x = x_mid;
                end

                if (y_mid * y_lower < 0)
                    x_upper = x_mid;
                else
                    x_lower = x_mid;

                iter_count = iter_count + 1;
            end
            if (findzeropoint_rect_x == -1)
                $display("No zero point found in the given range. Please try again with different range.");
        end
    endfunction

    analog function real findzeropoint_rect_y;
        input x_lower_input, x_upper_input, Vy, m;
        real x_lower_input, x_upper_input, Vy, m;
        real iter_count;
        real x_lower, x_upper;
        real y_mid;

        real y_lower, x_mid;

        begin
            x_lower = x_lower_input;
            x_upper = x_upper_input;
            iter_count = 0;
            y_mid = 1;
            findzeropoint_rect_y = -1;
            while (iter_count<`MAX_ITER && abs(y_mid) > `TOLERANCE) begin
                x_mid = (x_lower + x_upper) / 2;
                y_lower = tan(x_lower / 2 - m * `M_PI / 2) - sqrt((pow(Vy,2) - pow(x_lower,2))) / x_lower;
                y_mid = tan(x_mid / 2 - m * `M_PI / 2) - sqrt((pow(Vy,2) - pow(x_mid,2))) / x_mid;

                if (abs(y_mid) < `TOLERANCE)
                    findzeropoint_rect_y = x_mid;

                if (y_mid * y_lower < 0)
                    x_upper = x_mid;
                else
                    x_lower = x_mid;
                iter_count = iter_count + 1;
            end
            if (findzeropoint_rect_y == -1)
                $display("No zero point found in the given range. Please try again with different range.");
        end
    endfunction

    analog function real eim_rect;
        input lam, width, thickness, n_ps, n_sub, numflag, mode_index;
        real lam, width, thickness, n_ps, n_sub;
        integer numflag, mode_index;
        real k1, Vx, yx, hx, beta_x, neff_x_MMI, Vy, yy, hxy;
        real beta_xy[0:`MAX_MODES];
        real neff_xy[0:`MAX_MODES];
        real beta_xy_out[0:`MAX_MODES];
        real neff_xy_out[0:`MAX_MODES];
        integer i, m, mode_num;

        begin
            k1=2*`M_PI*n_ps/lam;
            Vx=2*`M_PI/lam*thickness*sqrt((pow(n_ps,2)-pow(n_sub,2)));
            yx=findzeropoint_rect_x(1e-6,min(`M_PI,Vx),Vx);
            hx=yx/thickness;
            beta_x=sqrt(pow(k1,2)-pow(hx,2));
            neff_x_MMI=lam*beta_x/(2*`M_PI);
            Vy=2*`M_PI/lam*width*sqrt((pow(neff_x_MMI,2)-pow(n_sub,2)));
            i=0;
            m=0;
            while (m*`M_PI<Vy) begin
                yy=findzeropoint_rect_y(m*`M_PI+1e-6,min((m+1)*`M_PI-1e-6,Vy),Vy,m);
                m=m+1;
                // if m*`M_PI>=Vy
                //     break
                hxy=yy/width;
                beta_xy[i]=sqrt(pow(beta_x,2)-pow(hxy,2));
                neff_xy[i]=lam*beta_xy[i]/(2*`M_PI);
                i=i+1;
            end
            mode_num=m-1;
            if (numflag==1 || numflag==2) begin
                for (i=0;i<mode_num;i=i+1) begin
                    beta_xy_out[i]=beta_xy[i];
                    neff_xy_out[i]=neff_xy[i];
                end
                if (numflag==1)
                    eim_rect=neff_xy_out[mode_index];
                else
                    eim_rect=beta_xy_out[mode_index];
            end
            else if (numflag==3) begin
                eim_rect=beta_x;
            end
            else if (numflag==4) begin
                eim_rect=mode_num;
            end
        end
    endfunction

    analog function real fieldfunc;
        input y,width,offset,ky,gamma,m,c0;
        real y,width,offset,ky,gamma,m,c0;
        real fai;

        begin
            fai = m * `M_PI / 2;
            if (((y - offset) <= (width / 2)) && ((-width / 2) <= (y - offset)))
                fieldfunc=c0 * cos(ky * (y - offset) - fai);
            else if ((y - offset) > width / 2)
                fieldfunc=c0 * cos(ky * width / 2 - fai) * exp(-gamma * ((y - offset) - width / 2));
            else
                fieldfunc=c0 * cos(ky * width / 2 + fai) * exp(gamma * ((y - offset) + width / 2));
        end
    endfunction

    analog function real normpower_fieldfunc_quad;
        input width,offset,ky,gamma,m,c0,width_MMI,border;
        real width,offset,ky,gamma,m,c0,width_MMI,border;
        real n,h;
        // real border=width_MMI/2;
        integer i;
        real sum;
        begin
            n=`IDX_NUM;
            h=border*2/n;
            sum=0.5*(pow(abs(fieldfunc(-1*border,width,offset,ky,gamma,m,c0)),2)+pow(abs(fieldfunc(border,width,offset,ky,gamma,m,c0)),2));
            for (i=1;i<n;i=i+1) begin
                sum=sum+pow(abs(fieldfunc(-1*border+i*h,width,offset,ky,gamma,m,c0)),2);
            end
            normpower_fieldfunc_quad=1/sqrt(h*sum);
        end
    endfunction

    analog function real infidx_modecoeffcient_quad;
        input width1,offset1,ky1,gamma1,m1,c01,width2,offset2,ky2,gamma2,m2,c02,border;
        real width1,offset1,ky1,gamma1,m1,c01,width2,offset2,ky2,gamma2,m2,c02,border;
        real n,h;
        // real border=20;
        integer i;
        real sum;
        begin
            n=`IDX_NUM;
            h=border*2/n;
            sum=0.5*(fieldfunc(-1*border,width1,offset1,ky1,gamma1,m1,c01)* fieldfunc(-1*border,width2,offset2,ky2,gamma2,m2,c02)+fieldfunc(border,width1,offset1,ky1,gamma1,m1,c01)* fieldfunc(border,width2,offset2,ky2,gamma2,m2,c02));
            for (i=1;i<n;i=i+1) begin
                sum=sum+fieldfunc(-1*border+i*h,width1,offset1,ky1,gamma1,m1,c01)* fieldfunc(-1*border+i*h,width2,offset2,ky2,gamma2,m2,c02);
            end
            infidx_modecoeffcient_quad=h*sum;
        end
    endfunction

    analog initial begin
        ulam0=design_lam*1e6;
        uwidth_in=width_port*1e6;
        uwidth_MMI=width_MMI*1e6;
        uoffset=port_offset*1e6;
        uthickness=thickness*1e6;
        alpha0=alpha_mmi*1e-6;
        k1_ref = 2 * `M_PI * n_ps / ulam0;
        kc_ref = 2 * `M_PI * n_sub / ulam0;

        mode_num_ref = eim_rect(ulam0, uwidth_MMI, uthickness, n_ps, n_sub, 4, 0);
        for (i=0;i<mode_num_ref;i=i+1) begin
            beta_MMI[i] = eim_rect(ulam0, uwidth_MMI, uthickness, n_ps, n_sub, 2, i);
        end
        length_MMI = p*1.5*`M_PI/(beta_MMI[0]-beta_MMI[1]); //1.5Lpi?
        // $display("length_MMI=%f",length_MMI);
        border=20;
    end

    analog begin
        ulam = V(Ilam1) * 1e6;
        beta_in=eim_rect(ulam, uwidth_in, uthickness, n_ps, n_sub, 2, 0);
        // $display("beta_in=%f",beta_in);
        beta_x_in=eim_rect(ulam, uwidth_in, uthickness, n_ps, n_sub, 3, 0);
        // $display("beta_x_in=%f",beta_x_in);
        beta_x_MMI=eim_rect(ulam, uwidth_MMI, uthickness, n_ps, n_sub, 3, 0);
        // $display("beta_x_MMI=%f",beta_x_MMI);
        mode_num = eim_rect(ulam, uwidth_MMI, uthickness, n_ps, n_sub, 4, 0);
        // $display("mode_num=%d",mode_num);
        ky_in = sqrt(pow(beta_x_in,2) - pow(beta_in,2));
        // $display("ky_in=%f",ky_in);
        gamma_in = sqrt(pow(beta_in,2) - pow(kc_ref,2));
        for (i=0;i<mode_num;i=i+1) begin
            neff_MMI[i] = eim_rect(ulam, uwidth_MMI, uthickness, n_ps, n_sub, 1, i);
            beta_MMI[i] = eim_rect(ulam, uwidth_MMI, uthickness, n_ps, n_sub, 2, i);
            ky_MMI[i] = sqrt(pow(beta_x_MMI,2) - pow(beta_MMI[i],2));
            gamma_MMI[i] = sqrt(pow(beta_MMI[i],2) - pow(kc_ref,2));
            // $display("neff_MMI[%d]=%f, beta_MMI[%d]=%f, ky_MMI[%d]=%f, gamma_MMI[%d]=%f",i,neff_MMI[i],i,beta_MMI[i],i,ky_MMI[i],i,gamma_MMI[i]);
        end

        c0=normpower_fieldfunc_quad(uwidth_in,0,ky_in,gamma_in,0,1,uwidth_MMI,border);
        // $display("c0=%f",c0);
        for (i=0;i<mode_num;i=i+1) begin
            c0_dis[i]=normpower_fieldfunc_quad(uwidth_MMI,0,ky_MMI[i],gamma_MMI[i],i,1,uwidth_MMI,border);
            // $display("c0_dis[%d]=%f",i,c0_dis[i]);
            icoff[i]=infidx_modecoeffcient_quad(uwidth_in,uoffset,ky_in,gamma_in,0,c0,uwidth_MMI,0,ky_MMI[i],gamma_MMI[i],i,c0_dis[i],border);
            // $display("icoff[%d]=%f",i,icoff[i]);
        end

        delta_beta = beta_MMI[0] - beta_MMI[1];
        // $display("delta_beta=%f",delta_beta);
        for (i=2;i<mode_num;i=i+1) begin
            beta_MMI[i]=beta_MMI[0] - i* (i + 2) * delta_beta / 3;
            // $display("beta_MMI[%d]=%f",i,beta_MMI[i]);
        end
        for (i=0;i<mode_num;i=i+1) begin
            ocoff_r[i]=icoff[i]*cos(-beta_MMI[i]*length_MMI);
            ocoff_i[i]=icoff[i]*sin(-beta_MMI[i]*length_MMI)*-1.0;
            // $display("ocoff_r[%d]=%f, ocoff_i[%d]=%f",i,ocoff_r[i],i,ocoff_i[i]);
        end
        h=border*2.0/`IDX_NUM;
        sum_outfield1=0;
        sum_outfield2=0;
        sum_outfield3=0;
        sum_outfield4=0;
        for (j=0;j<mode_num;j=j+1) begin
            sum_singlemode1=0.5*(ocoff_r[j]*fieldfunc(-1*border,uwidth_MMI,0,ky_MMI[j],gamma_MMI[j],j,c0_dis[j])*fieldfunc(-1*border,uwidth_in,uoffset,ky_in,gamma_in,0,c0)+ocoff_r[j]*fieldfunc(border,uwidth_MMI,0,ky_MMI[j],gamma_MMI[j],j,c0_dis[j])*fieldfunc(border,uwidth_in,uoffset,ky_in,gamma_in,0,c0));
            sum_singlemode2=0.5*(ocoff_i[j]*fieldfunc(-1*border,uwidth_MMI,0,ky_MMI[j],gamma_MMI[j],j,c0_dis[j])*fieldfunc(-1*border,uwidth_in,uoffset,ky_in,gamma_in,0,c0)+ocoff_i[j]*fieldfunc(border,uwidth_MMI,0,ky_MMI[j],gamma_MMI[j],j,c0_dis[j])*fieldfunc(border,uwidth_in,uoffset,ky_in,gamma_in,0,c0));
            sum_singlemode3=0.5*(ocoff_r[j]*fieldfunc(-1*border,uwidth_MMI,0,ky_MMI[j],gamma_MMI[j],j,c0_dis[j])*fieldfunc(-1*border,uwidth_in,-uoffset,ky_in,gamma_in,0,c0)+ocoff_r[j]*fieldfunc(border,uwidth_MMI,0,ky_MMI[j],gamma_MMI[j],j,c0_dis[j])*fieldfunc(border,uwidth_in,-uoffset,ky_in,gamma_in,0,c0));
            sum_singlemode4=0.5*(ocoff_i[j]*fieldfunc(-1*border,uwidth_MMI,0,ky_MMI[j],gamma_MMI[j],j,c0_dis[j])*fieldfunc(-1*border,uwidth_in,-uoffset,ky_in,gamma_in,0,c0)+ocoff_i[j]*fieldfunc(border,uwidth_MMI,0,ky_MMI[j],gamma_MMI[j],j,c0_dis[j])*fieldfunc(border,uwidth_in,-uoffset,ky_in,gamma_in,0,c0));
            for (i=1;i<`IDX_NUM;i=i+1) begin
                sum_singlemode1=sum_singlemode1+ocoff_r[j]*fieldfunc(-1*border+i*h,uwidth_MMI,0,ky_MMI[j],gamma_MMI[j],j,c0_dis[j])*fieldfunc(-1*border+i*h,uwidth_in,uoffset,ky_in,gamma_in,0,c0);
                sum_singlemode2=sum_singlemode2+ocoff_i[j]*fieldfunc(-1*border+i*h,uwidth_MMI,0,ky_MMI[j],gamma_MMI[j],j,c0_dis[j])*fieldfunc(-1*border+i*h,uwidth_in,uoffset,ky_in,gamma_in,0,c0);
                sum_singlemode3=sum_singlemode3+ocoff_r[j]*fieldfunc(-1*border+i*h,uwidth_MMI,0,ky_MMI[j],gamma_MMI[j],j,c0_dis[j])*fieldfunc(-1*border+i*h,uwidth_in,-uoffset,ky_in,gamma_in,0,c0);
                sum_singlemode4=sum_singlemode4+ocoff_i[j]*fieldfunc(-1*border+i*h,uwidth_MMI,0,ky_MMI[j],gamma_MMI[j],j,c0_dis[j])*fieldfunc(-1*border+i*h,uwidth_in,-uoffset,ky_in,gamma_in,0,c0);
            end
            // $display("sum_singlemode1=%f, sum_singlemode2=%f, sum_singlemode3=%f, sum_singlemode4=%f",sum_singlemode1,sum_singlemode2,sum_singlemode3,sum_singlemode4);
            sum_outfield1=sum_outfield1+h*sum_singlemode1;
            sum_outfield2=sum_outfield2+h*sum_singlemode2;
            sum_outfield3=sum_outfield3+h*sum_singlemode3;
            sum_outfield4=sum_outfield4+h*sum_singlemode4;
        end
        c1r=sum_outfield1;
        c1i=sum_outfield2;
        c2r=sum_outfield3;
        c2i=sum_outfield4;
        // $display("c1r=%f, c1i=%f, c2r=%f, c2i=%f",c1r,c1i,c2r,c2i);

        in1r=sqrt(V(Ipow1))*cos(V(Iphase1)/360.0*2*`M_PI);
        in1i=sqrt(V(Ipow1))*sin(V(Iphase1)/360.0*2*`M_PI);
        in2r=sqrt(V(Ipow2))*cos(V(Iphase2)/360.0*2*`M_PI);
        in2i=sqrt(V(Ipow2))*sin(V(Iphase2)/360.0*2*`M_PI);
        out1r=c1r*in1r-c1i*in1i+c2r*in2r-c2i*in2i;
        out1i=c1r*in1i+c1i*in1r+c2r*in2i+c2i*in2r;
        out2r=c2r*in1r-c2i*in1i-c1r*in2r+c1i*in2i;
        out2i=c2r*in1i+c2i*in1r-c1r*in2i-c1i*in2r;

        V(Opow1)   <+ pow(out1r,2)+pow(out1i,2);
        V(Ophase1) <+ atan2(out1i,out1r)*360.0/(2*`M_PI);
        V(Olam1)   <+ ulam/1e6;
        V(Opow2)   <+ pow(out2r,2)+pow(out2i,2);
        V(Ophase2) <+ atan2(out2i,out2r)*360.0/(2*`M_PI);
        V(Olam2)   <+ ulam/1e6;
    end

endmodule