//VerilogA for modulator,basicea,veriloga
`include "constants.vams"
`include "disciplines.vams"

`define MAX_ARRAY 1000

module basicea(V1,V2,Gnd,Ipow,Iphase,Ilam,Opow,Ophase,Olam);
	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 Ipow,Iphase,Ilam;
	output Opow,Ophase,Olam;
	inout V1,V2,Gnd;
	electrical V1,V2,Ipow,Iphase,Ilam,Opow,Ophase,Olam,Gnd,Nph,Ns;
	
	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];

	branch (Gnd,Ns) Cs,Is;
	branch (Ns,Gnd) It;
	branch (Gnd,Nph) Cph,Iph;

	analog begin
		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(V1)-V(V2);
		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),i); //pow(N,i)*pow(f,i)
				Alpha=Alpha+alpha_array[i*num_carriler+j]*pow(Vcalc,j)*pow(V(Ns),i); //Caution the unit
			end
		end
		if (Gamma==0)
			opow=0;
		else begin
			opow=pow(sqrt(V(Ipow))*exp(-Gamma/2),2);
		end
		
		I(Is) <+ f*V(Ilam)/(`P_H*`P_C)*exp(-Gam_c)*(1-exp(-Gamma+2*Gam_c))*V(Ipow);
		I(It) <+ V(Ns)/tau;
		I(Cs) <+ ddt(V(Cs));

		I(Iph) <+ Alpha/2*ddt(Gamma);
		I(Cph) <+ ddt(V(Cph));
		ophase=V(Nph)+V(Iphase)*(2*`M_PI)/360.0;

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

endmodule