//VerilogA for passive,xcouplerbasic,veriloga

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

module xcouplerbasic(Ipow1,Iphase1,Ilam1,Ipow2,Iphase2,Ilam2,Opow1,Ophase1,Olam1,Opow2,Ophase2,Olam2);
    parameter integer modespec = 1 from [1:2];
    parameter integer alpha_use_dBm = 1 from [0:1];
    parameter real reffreq = 193.1e12 from (0:inf);
    parameter real l_upper = 1e-3 from [0:inf);
    parameter real l_lower = 1e-3 from [0:inf);
    parameter real l_coupler = 392.7e-6 from [0:inf);
    parameter real neff_upper_te = 2.6 from (0:inf);
    parameter real neff_upper_tm = 2.6 from (0:inf);
    parameter real neff_lower_te = 2.6 from (0:inf);
    parameter real neff_lower_tm = 2.6 from (0:inf);
    parameter real n_gv_upper_te = 4.2 from (0:inf);
    parameter real n_gv_upper_tm = 4.2 from (0:inf);
    parameter real n_gv_lower_te = 4.2 from (0:inf);
    parameter real n_gv_lower_tm = 4.2 from (0:inf);
    parameter real disper_upper_te = 0 from (-inf:inf);
    parameter real disper_upper_tm = 0 from (-inf:inf);
    parameter real disper_lower_te = 0 from (-inf:inf);
    parameter real disper_lower_tm = 0 from (-inf:inf);
    parameter real attenu_upper_te = 0 from [0:inf);
    parameter real attenu_upper_tm = 0 from [0:inf);
    parameter real attenu_lower_te = 0 from [0:inf);
    parameter real attenu_lower_tm = 0 from [0:inf);
    parameter real coupling_te = 2000 from (-inf:inf);
    parameter real coupling_tm = 2000 from (-inf: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 eeir1,eeii1,emir1,emii1;
    real eeir2,eeii2,emir2,emii2;
    real eeor1,eeoi1,emor1,emoi1;
    real eeor2,eeoi2,emor2,emoi2;

    real r_sqrt_te,i_sqrt_te,a_sqrt_te,t_sqrt_te;
    real r_sqrt_tm,i_sqrt_tm,a_sqrt_tm,t_sqrt_tm;
    real r_clow;
    real z1r_te,z1i_te,z21r_te,z21i_te,z2r_te,z2i_te,z3r_te,z3i_te;
    real z1r_tm,z1i_tm,z21r_tm,z21i_tm,z2r_tm,z2i_tm,z3r_tm,z3i_tm;
    real c11r_te,c11i_te,c12r_te,c12i_te,c21r_te,c21i_te,c22r_te,c22i_te,c31r_te,c31i_te,c32r_te,c32i_te;
    real c11r_tm,c11i_tm,c12r_tm,c12i_tm,c21r_tm,c21i_tm,c22r_tm,c22i_tm,c31r_tm,c31i_tm,c32r_tm,c32i_tm;
    real c1r_te,c1i_te,c2r_te,c2i_te,c3r_te,c3i_te;
    real c1r_tm,c1i_tm,c2r_tm,c2i_tm,c3r_tm,c3i_tm;
    real t11r,t11i,t12r,t12i,t13r,t13i,t14r,t14i,t21r,t21i,t22r,t22i,t23r,t23i,t24r,t24i;
    real t31r,t31i,t32r,t32i,t33r,t33i,t34r,t34i,t41r,t41i,t42r,t42i,t43r,t43i,t44r,t44i;

    real alpha_upper_te,alpha_upper_tm,alpha_lower_te,alpha_lower_tm;
    real beta_upper_ter,beta_upper_tei,beta_upper_tmr,beta_upper_tmi,beta_lower_ter,beta_lower_tei,beta_lower_tmr,beta_lower_tmi;
    real beta_avg_ter,beta_avg_tei,beta_avg_tmr,beta_avg_tmi;
    real beta_diff_ter,beta_diff_tei,beta_diff_tmr,beta_diff_tmi;
    real delta_ter,delta_tei,delta_tmr,delta_tmi;
    real theta_ter,theta_tei,theta_tmr,theta_tmi;
    real r_ter,r_tei,r_tmr,r_tmi;
    real s_ter,s_tei,s_tmr,s_tmi;
    real freq_upper,freq_lower;
    real iph_upper,iph_lower;

    analog begin
        iph_upper=V(Iphase1)/360.0*2*`M_PI;
        iph_lower=V(Iphase2)/360.0*2*`M_PI;
        if (V(Ilam1)==0)
            freq_upper=`P_C/V(Ilam2);
        else
            freq_upper=`P_C/V(Ilam1);
        if (V(Ilam2)==0)
            freq_lower=`P_C/V(Ilam1);
        else
            freq_lower=`P_C/V(Ilam2);
        if (modespec ==1) begin
            eeir1=sqrt(V(Ipow1))*cos(iph_upper);
            eeii1=sqrt(V(Ipow1))*sin(iph_upper);
            emir1=0;
            emii1=0;
            eeir2=sqrt(V(Ipow2))*cos(iph_lower);
            eeii2=sqrt(V(Ipow2))*sin(iph_lower);
            emir2=0;
            emii2=0;
        end
        else if (modespec ==2) begin
            eeir1=0;
            eeii1=0;
            emir1=sqrt(V(Ipow1))*cos(iph_upper);
            emii1=sqrt(V(Ipow1))*sin(iph_upper);
            eeir2=0;
            eeii2=0;
            emir2=sqrt(V(Ipow2))*cos(iph_lower);
            emii2=sqrt(V(Ipow2))*sin(iph_lower);
        end

        if (alpha_use_dBm == 1) begin
            alpha_upper_te=0.1*ln(10)*attenu_upper_te;
            alpha_upper_tm=0.1*ln(10)*attenu_upper_tm;
            alpha_lower_te=0.1*ln(10)*attenu_lower_te;
            alpha_lower_tm=0.1*ln(10)*attenu_lower_tm;
        end
        else begin
            alpha_upper_te=attenu_upper_te;
            alpha_upper_tm=attenu_upper_tm;
            alpha_lower_te=attenu_lower_te;
            alpha_lower_tm=attenu_lower_tm;
        end

        beta_upper_ter=2.0*`M_PI*reffreq/`P_C*neff_upper_te+2*`M_PI/`P_C*n_gv_upper_te*(freq_upper-reffreq)-`M_PI*`P_C*disper_upper_te/pow(reffreq,2)*pow(freq_upper-reffreq,2);
        beta_upper_tei=-0.5*alpha_upper_te;
        beta_upper_tmr=2.0*`M_PI*reffreq/`P_C*neff_upper_tm+2*`M_PI/`P_C*n_gv_upper_tm*(freq_upper-reffreq)-`M_PI*`P_C*disper_upper_tm/pow(reffreq,2)*pow(freq_upper-reffreq,2);
        beta_upper_tmi=-0.5*alpha_upper_tm;
        beta_lower_ter=2.0*`M_PI*reffreq/`P_C*neff_lower_te+2*`M_PI/`P_C*n_gv_lower_te*(freq_lower-reffreq)-`M_PI*`P_C*disper_lower_te/pow(reffreq,2)*pow(freq_lower-reffreq,2);
        beta_lower_tei=-0.5*alpha_lower_te;
        beta_lower_tmr=2.0*`M_PI*reffreq/`P_C*neff_lower_tm+2*`M_PI/`P_C*n_gv_lower_tm*(freq_lower-reffreq)-`M_PI*`P_C*disper_lower_tm/pow(reffreq,2)*pow(freq_lower-reffreq,2);
        beta_lower_tmi=-0.5*alpha_lower_tm;

        beta_avg_ter=(beta_upper_ter+beta_lower_ter)/2.0;
        beta_avg_tei=(beta_upper_tei+beta_lower_tei)/2.0;
        beta_avg_tmr=(beta_upper_tmr+beta_lower_tmr)/2.0;
        beta_avg_tmi=(beta_upper_tmi+beta_lower_tmi)/2.0;
        beta_diff_ter=(beta_upper_ter-beta_lower_ter)/2.0;
        beta_diff_tei=(beta_upper_tei-beta_lower_tei)/2.0;
        beta_diff_tmr=(beta_upper_tmr-beta_lower_tmr)/2.0;
        beta_diff_tmi=(beta_upper_tmi-beta_lower_tmi)/2.0;

        r_sqrt_te=pow(beta_diff_ter,2)+pow(beta_diff_tei,2)+pow(coupling_te,2);
        i_sqrt_te=2.0*beta_diff_ter*beta_diff_tei;
        a_sqrt_te=sqrt(pow(r_sqrt_te,2)+pow(i_sqrt_te,2));
        t_sqrt_te=atan2(i_sqrt_te,r_sqrt_te);
        delta_ter=sqrt(a_sqrt_te)*cos(t_sqrt_te/2);
        delta_tei=sqrt(a_sqrt_te)*sin(t_sqrt_te/2);
        r_sqrt_tm=pow(beta_diff_tmr,2)+pow(beta_diff_tmi,2)+pow(coupling_tm,2);
        i_sqrt_tm=2.0*beta_diff_tmr*beta_diff_tmi;
        a_sqrt_tm=sqrt(pow(r_sqrt_tm,2)+pow(i_sqrt_tm,2));
        t_sqrt_tm=atan2(i_sqrt_tm,r_sqrt_tm);
        delta_tmr=sqrt(a_sqrt_tm)*cos(t_sqrt_tm/2);
        delta_tmi=sqrt(a_sqrt_tm)*sin(t_sqrt_tm/2);

        theta_ter=delta_ter*l_coupler;
        theta_tei=delta_tei*l_coupler;
        theta_tmr=delta_tmr*l_coupler;
        theta_tmi=delta_tmi*l_coupler;

        r_clow=pow(delta_ter,2)+pow(delta_tei,2);
        r_ter=(beta_diff_ter*delta_ter+beta_diff_tei*delta_tei)/r_clow;
        r_tei=(-beta_diff_ter*delta_tei+beta_diff_tei*delta_ter)/r_clow;
        r_tmr=(beta_diff_tmr*delta_tmr+beta_diff_tmi*delta_tmi)/r_clow;
        r_tmi=(-beta_diff_tmr*delta_tmi+beta_diff_tmi*delta_tmr)/r_clow;

        s_ter=coupling_te*delta_ter/(pow(delta_ter,2)+pow(delta_tei,2));
        s_tei=-coupling_te*delta_tei/(pow(delta_ter,2)+pow(delta_tei,2));
        s_tmr=coupling_tm*delta_tmr/(pow(delta_tmr,2)+pow(delta_tmi,2));
        s_tmi=-coupling_tm*delta_tmi/(pow(delta_tmr,2)+pow(delta_tmi,2));

        //cos(theta_te)
        z1r_te=(exp(-theta_tei)+exp(theta_tei))/2.0*cos(theta_ter);
        z1i_te=(exp(-theta_tei)-exp(theta_tei))/2.0*sin(theta_ter);
        //sin(theta_te)
        z21r_te=(exp(theta_tei)+exp(-theta_tei))/2.0*sin(theta_ter);
        z21i_te=(exp(theta_tei)-exp(-theta_tei))/2.0*cos(theta_ter);
        //R*sin(theta_te)
        z2r_te=r_ter*z21r_te-r_tei*z21i_te;
        z2i_te=r_ter*z21i_te+r_tei*z21r_te;
        //S*sin(theta_te)
        z3r_te=s_ter*z21r_te-s_tei*z21i_te;
        z3i_te=s_ter*z21i_te+s_tei*z21r_te;
        //cos(theta_tm)
        z1r_tm=(exp(-theta_tmi)+exp(theta_tmi))/2.0*cos(theta_tmr);
        z1i_tm=(exp(-theta_tmi)-exp(theta_tmi))/2.0*sin(theta_tmr);
        //sin(theta_tm)
        z21r_tm=(exp(theta_tmi)+exp(-theta_tmi))/2.0*sin(theta_tmr);
        z21i_tm=(exp(theta_tmi)-exp(-theta_tmi))/2.0*cos(theta_tmr);
        //R*sin(theta_tm)
        z2r_tm=r_tmr*z21r_tm-r_tmi*z21i_tm;
        z2i_tm=r_tmr*z21i_tm+r_tmi*z21r_tm;
        //S*sin(theta_tm)
        z3r_tm=s_tmr*z21r_tm-s_tmi*z21i_tm;
        z3i_tm=s_tmr*z21i_tm+s_tmi*z21r_tm;

        c11r_te=z1r_te+z2i_te;
        c11i_te=z1i_te-z2r_te;
        c12r_te=exp(beta_avg_tei*l_coupler+beta_upper_tei*(l_upper-l_coupler))*cos(-beta_avg_ter*l_coupler-beta_upper_ter*(l_upper-l_coupler));
        c12i_te=exp(beta_avg_tei*l_coupler+beta_upper_tei*(l_upper-l_coupler))*sin(-beta_avg_ter*l_coupler-beta_upper_ter*(l_upper-l_coupler));
        c21r_te=z1r_te-z2i_te;
        c21i_te=z1i_te+z2r_te;
        c22r_te=exp(beta_avg_tei*l_coupler+beta_lower_tei*(l_lower-l_coupler))*cos(-beta_avg_ter*l_coupler-beta_lower_ter*(l_lower-l_coupler));
        c22i_te=exp(beta_avg_tei*l_coupler+beta_lower_tei*(l_lower-l_coupler))*sin(-beta_avg_ter*l_coupler-beta_lower_ter*(l_lower-l_coupler));
        c31r_te=z3i_te;
        c31i_te=-z3r_te;
        c32r_te=exp(0.5*(beta_upper_tei*l_upper+beta_lower_tei*l_lower))*cos(-0.5*(beta_upper_ter*l_upper+beta_lower_ter*l_lower));
        c32i_te=exp(0.5*(beta_upper_tei*l_upper+beta_lower_tei*l_lower))*sin(-0.5*(beta_upper_ter*l_upper+beta_lower_ter*l_lower));
        c11r_tm=z1r_tm+z2i_tm;
        c11i_tm=z1i_tm-z2r_tm;
        c12r_tm=exp(beta_avg_tmi*l_coupler+beta_upper_tmi*(l_upper-l_coupler))*cos(-beta_avg_tmr*l_coupler-beta_upper_tmr*(l_upper-l_coupler));
        c12i_tm=exp(beta_avg_tmi*l_coupler+beta_upper_tmi*(l_upper-l_coupler))*sin(-beta_avg_tmr*l_coupler-beta_upper_tmr*(l_upper-l_coupler));
        c21r_tm=z1r_tm-z2i_tm;
        c21i_tm=z1i_tm+z2r_tm;
        c22r_tm=exp(beta_avg_tmi*l_coupler+beta_lower_tmi*(l_lower-l_coupler))*cos(-beta_avg_tmr*l_coupler-beta_lower_tmr*(l_lower-l_coupler));
        c22i_tm=exp(beta_avg_tmi*l_coupler+beta_lower_tmi*(l_lower-l_coupler))*sin(-beta_avg_tmr*l_coupler-beta_lower_tmr*(l_lower-l_coupler));
        c31r_tm=z3i_tm;
        c31i_tm=-z3r_tm;
        c32r_tm=exp(0.5*(beta_upper_tmi*l_upper+beta_lower_tmi*l_lower))*cos(-0.5*(beta_upper_tmr*l_upper+beta_lower_tmr*l_lower));
        c32i_tm=exp(0.5*(beta_upper_tmi*l_upper+beta_lower_tmi*l_lower))*sin(-0.5*(beta_upper_tmr*l_upper+beta_lower_tmr*l_lower));

        c1r_te=c11r_te*c12r_te-c11i_te*c12i_te;
        c1i_te=c11r_te*c12i_te+c11i_te*c12r_te;
        c2r_te=c21r_te*c22r_te-c21i_te*c22i_te;
        c2i_te=c21r_te*c22i_te+c21i_te*c22r_te;
        c3r_te=c31r_te*c32r_te-c31i_te*c32i_te;
        c3i_te=c31r_te*c32i_te+c31i_te*c32r_te;
        c1r_tm=c11r_tm*c12r_tm-c11i_tm*c12i_tm;
        c1i_tm=c11r_tm*c12i_tm+c11i_tm*c12r_tm;
        c2r_tm=c21r_tm*c22r_tm-c21i_tm*c22i_tm;
        c2i_tm=c21r_tm*c22i_tm+c21i_tm*c22r_tm;
        c3r_tm=c31r_tm*c32r_tm-c31i_tm*c32i_tm;
        c3i_tm=c31r_tm*c32i_tm+c31i_tm*c32r_tm;

        t11r=c1r_te;
        t11i=c1i_te;
        t12r=0;
        t12i=0;
        t13r=c3r_te;
        t13i=c3i_te;
        t14r=0;
        t14i=0;
        t21r=0;
        t21i=0;
        t22r=c1r_tm;
        t22i=c1i_tm;
        t23r=0;
        t23i=0;
        t24r=c3r_tm;
        t24i=c3i_tm;
        t31r=c3r_te;
        t31i=c3i_te;
        t32r=0;
        t32i=0;
        t33r=c2r_te;
        t33i=c2i_te;
        t34r=0;
        t34i=0;
        t41r=0;
        t41i=0;
        t42r=c3r_tm;
        t42i=c3i_tm;
        t43r=0;
        t43i=0;
        t44r=c2r_tm;
        t44i=c2i_tm;

        eeor1=t11r*eeir1-t11i*eeii1+t12r*emir1-t12i*emii1+t13r*eeir2-t13i*eeii2+t14r*emir2-t14i*emii2;
        eeoi1=t11r*eeii1+t11i*eeir1+t12r*emii1+t12i*emir1+t13r*eeii2+t13i*eeir2+t14r*emii2+t14i*emir2;
        emor1=t21r*eeir1-t21i*eeii1+t22r*emir1-t22i*emii1+t23r*eeir2-t23i*eeii2+t24r*emir2-t24i*emii2;
        emoi1=t21r*eeii1+t21i*eeir1+t22r*emii1+t22i*emir1+t23r*eeii2+t23i*eeir2+t24r*emii2+t24i*emir2;
        eeor2=t31r*eeir1-t31i*eeii1+t32r*emir1-t32i*emii1+t33r*eeir2-t33i*eeii2+t34r*emir2-t34i*emii2;
        eeoi2=t31r*eeii1+t31i*eeir1+t32r*emii1+t32i*emir1+t33r*eeii2+t33i*eeir2+t34r*emii2+t34i*emir2;
        emor2=t41r*eeir1-t41i*eeii1+t42r*emir1-t42i*emii1+t43r*eeir2-t43i*eeii2+t44r*emir2-t44i*emii2;
        emoi2=t41r*eeii1+t41i*eeir1+t42r*emii1+t42i*emir1+t43r*eeii2+t43i*eeir2+t44r*emii2+t44i*emir2;

        if (V(Ilam1)==0)
            V(Olam1) <+ V(Ilam2);
        else
            V(Olam1) <+ V(Ilam1);
        if (V(Ilam2)==0)
            V(Olam2) <+ V(Ilam1); 
        else
            V(Olam2) <+ V(Ilam2); 
        if (modespec == 1) begin
            V(Opow1) <+ pow(eeor1,2)+pow(eeoi1,2);
            V(Ophase1) <+ atan2(eeoi1,eeor1)*360.0/(2*`M_PI);
            V(Opow2) <+ pow(eeor2,2)+pow(eeoi2,2);
            V(Ophase2) <+ atan2(eeoi2,eeor2)*360.0/(2*`M_PI);
        end
        else if (modespec == 2) begin
            V(Opow1) <+ pow(emor1,2)+pow(emoi1,2);
            V(Ophase1) <+ atan2(emoi1,emor1)*360.0/(2*`M_PI);
            V(Opow2) <+ pow(emor2,2)+pow(emoi2,2);
            V(Ophase2) <+ atan2(emoi2,emor2)*360.0/(2*`M_PI);
        end
        
    end

endmodule