clear
clc

function todBconv = todB (v)
  todBconv = 10 * log10(v);
endfunction

function ser = Pse(x)
ser = 2 * qfunc(sqrt(x)) - qfunc(sqrt(x))^2;
endfunction

LightSpeed = 3e8;
k = - 228.6;
Tm = 275;

Fc = 12e9;
lambda = LightSpeed / Fc;

%%%%%% Satellite %%%%%%%
P = 5;
Lt = 1.25;
Theta3todBt = 2.8;
etat = 0.67;

Dt = 70 * lambda / Theta3todBt;
Gtmax = etat * (pi * Dt / lambda)^2;
Gt = Gtmax / 2;

disp("\nEIRP (dB)\n");
EIRP = todB(P) - todB(Lt) + todB(Gt)

%%%%%%% Medium %%%%%%%%%%%
d = 38e6;
La = 0.3;
Tm = 275;
Lrain = 1;

Lfs = (4 * pi * d / lambda)^2;

Tawr = Tm * (1 - 1 / 10^(La/10));

disp("\nFree space loss (dB)\n");
todB(Lfs)

disp("\nMedium losses without rain (dB)\n");
Lwor = todB(Lfs) + La

Tar = Tm * (1 - 1 / 10^( (La+Lrain)/10) );

disp("\nMedium losses with rain (dB)\n");
Lwr = Lwor + Lrain

DTa = Tar - Tawr;

%%%%%% Earth Station %%%%%%%
Dr = 2.2;
etar = 0.54;
Ts = 200;
Bn = 36e6;
Rb = 70e6;
alpha = 0.4;

Bw = (1+alpha) * Bn;
disp("\nNyquist filter bandwidth (MHz)\n");
Bw /1e6

Rs= Bn;
disp("\nSymbol rate (MSymbols/sec)\n");
Rs/1e6

disp("\nModulation order\n");
M = pow2(ceil(log2(2^(Rb/Rs))))

Grmax = etar * (pi * Dr / lambda)^2;

disp("\nG/T without rain (dB1/K)\n");
GoverT = todB(Grmax) - todB(Ts)

%%%% Results without rain %%%%
disp("\nC/N0 without rain (dBHz)\n");
CoverN0 = -k + EIRP - Lwor + GoverT

disp("\nC/N without rain (dB)\n");
CNRwr = CoverN0 - todB(Bn)

Psewr = Pse( 10^(CNRwr/10) );

disp("\nQPSK bit error probability (without rain)\n");
BERwr = Psewr / 2

%%%% Results with rain %%%%
GoverTr = todB(Grmax) - todB(Ts + DTa);

disp("\nC/N0 with rain (dBHz)\n");
CoverN0r = -k + EIRP - Lwr + GoverTr

disp("\nC/N with rain (dB)\n");
CNRr = CoverN0r - todB(Bn)

Pser = Pse( 10^(CNRr/10) );

disp("\nQPSK bit error probability (with rain)\n");
BERr = Pser / 2