//VerilogA for passive,wgstrmodecp,veriloga

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

module wgstrmodecp(Ipow,Iphase,Ilam,Opow,Ophase,Olam);
    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 length = 1e-3 from [0:inf);
    parameter real neff_te = 2.6 from (0:inf);
    parameter real neff_tm = 2.6 from (0:inf);
    parameter real n_gv_te = 4.2 from (0:inf);
    parameter real n_gv_tm = 4.2 from (0:inf);
    parameter real disper_te = 0;
    parameter real disper_tm = 0;
    parameter real attenu_te = 0 from [0:inf);//attenuation in dB/m
    parameter real attenu_tm = 0 from [0:inf);//attenuation in dB/m
    parameter real phasedelay_te = 0;
    parameter real phasedelay_tm = 0;
    parameter real coupling_tem = 0;

    input Ipow,Iphase,Ilam;
    output Opow,Ophase,Olam;
    electrical Ipow,Iphase,Ilam,Opow,Ophase,Olam;

    real eeir,eeii,emir,emii;
    real eeor,eeoi,emor,emoi;
    real t11r,t11i,t12r,t12i,t21r,t21i,t22r,t22i;
    real beta_ter,beta_tei,beta_tmr,beta_tmi,beta_avr,beta_avi,beta_dvr,beta_dvi;
    real rfront_te,rfront_tm,t11r1,t11i1,t1i2,t22r1,t22i1;
    real c1r,c1i,a,t,deltar,deltai,thetar,thetai,sr,si,rr,ri,z1r,z1i,z21r,z21i,z2r,z2i,z3r,z3i;
    real c11r,c11i,c12r,c12i,c21r,c21i,c22r,c22i,c31r,c31i;
    real phasedelay_te_u=phasedelay_te/360.0*2*`M_PI;
    real phasedelay_tm_u=phasedelay_tm/360.0*2*`M_PI;
    real alpha_te,alpha_tm;
    real freq;
    real iph;

    analog begin
        freq=`P_C/V(Ilam);
        iph=V(Iphase)/360.0*2*`M_PI;
        if (modespec == 1) begin
            eeir=sqrt(V(Ipow))*cos(iph);
            eeii=sqrt(V(Ipow))*sin(iph);
            emir=0;
            emii=0;
        end
        else if (modespec == 2) begin
            eeir=0;
            eeii=0;
            emir=sqrt(V(Ipow))*cos(iph);
            emii=sqrt(V(Ipow))*sin(iph);
        end

        if (alpha_use_dBm == 1) begin
            alpha_te=0.1*ln(10)*attenu_te;
            alpha_tm=0.1*ln(10)*attenu_tm;
        end
        else begin
            alpha_te=attenu_te;
            alpha_tm=attenu_tm;
        end

        beta_ter=2*`M_PI*reffreq/`P_C*neff_te+2*`M_PI/`P_C*n_gv_te*(freq-reffreq)-`M_PI*`P_C*disper_te/pow(reffreq,2)*pow(freq-reffreq,2);     
        beta_tei=-alpha_te/2;
        beta_tmr=2*`M_PI*reffreq/`P_C*neff_tm+2*`M_PI/`P_C*n_gv_tm*(freq-reffreq)-`M_PI*`P_C*disper_tm/pow(reffreq,2)*pow(freq-reffreq,2);
        beta_tmi=-alpha_tm/2;

        beta_avr=(beta_ter+beta_tmr)/2;
        beta_avi=(beta_tmi+beta_tmi)/2;
        beta_dvr=(beta_ter-beta_tmr)/2;
        beta_dvi=(beta_tmi-beta_tmi)/2;
        c1r=beta_dvr*beta_dvr+beta_dvi*beta_dvi+coupling_tem*coupling_tem;
        c1i=2*beta_dvr*beta_dvi;
        a=sqrt(pow(c1r,2)+pow(c1i,2));
        t=atan2(c1i,c1r);
        deltar=sqrt(a)*cos(t/2);
        deltai=sqrt(a)*sin(t/2);
        thetar=deltar*length;
        thetai=deltai*length;
        sr=coupling_tem*deltar/(deltar*deltar+deltai*deltai);
        si=-coupling_tem*deltai/(deltar*deltar+deltai*deltai);
        rr=(beta_dvr*deltar+beta_dvi*deltai)/(deltar*deltar+deltai*deltai);
        ri=(-beta_dvr*deltai+beta_dvi*deltar)/(deltar*deltar+deltai*deltai);

        //cos(theta)
        z1r=(exp(-thetai)+exp(thetai))/2.0*cos(thetar);
        z1i=(exp(-thetai)-exp(thetai))/2.0*sin(thetar);
        //sin(theta)
        z21r=(exp(thetai)+exp(-thetai))/2.0*sin(thetar);
        z21i=(exp(thetai)-exp(-thetai))/2.0*cos(thetar);
        //R*sin(theta)
        z2r=rr*z21r-ri*z21i;
        z2i=rr*z21i+ri*z21r;
        //S*sin(theta)
        z3r=sr*z21r-si*z21i;
        z3i=sr*z21i+si*z21r;

        c11r=z1r+z2i;
        c11i=z1i-z2r;
        c12r=exp(beta_avi*length)*cos(-phasedelay_te_u-beta_avr*length);
        c12i=exp(beta_avi*length)*sin(-phasedelay_te_u-beta_avr*length);
        c21r=z1r-z2i;
        c21i=z1i+z2r;
        c22r=exp(beta_avi*length)*cos(-phasedelay_tm_u-beta_avr*length);
        c22i=exp(beta_avi*length)*sin(-phasedelay_tm_u-beta_avr*length);
        c31r=z3i;
        c31i=-z3r;

        t11r=c11r*c12r-c11i*c12i;
        t11i=c11r*c12i+c11i*c12r;
        t12r=c31r*c12r-c31i*c12i;
        t12i=c31r*c12i+c31i*c12r;
        t21r=c31r*c22r-c31i*c22i;
        t21i=c31r*c22i+c31i*c22r;
        t22r=c21r*c22r-c21i*c22i;
        t22i=c21r*c22i+c21i*c22r;

        eeor=t11r*eeir-t11i*eeii+t12r*emir-t12i*emii;
        eeoi=t11r*eeii+t11i*eeir+t12r*emii+t12i*emir;
        emor=t21r*eeir-t21i*eeii+t22r*emir-t22i*emii;
        emoi=t21r*eeii+t21i*eeir+t22r*emii+t22i*emir;

        V(Olam) <+ V(Ilam);
        if (modespec == 1) begin
            V(Opow) <+ eeor*eeor+eeoi*eeoi;
            V(Ophase) <+ atan2(eeoi,eeor)*360.0/(2*`M_PI);
        end
        else if (modespec == 2) begin
            V(Opow) <+ emor*emor+emoi*emoi;
            V(Ophase) <+ atan2(emoi,emor)*360.0/(2*`M_PI);
        end
    end

endmodule