#
# mimo.r
#

# Konstanten und Parameter setzen
# Empfangsantennen (nicht ändern wegen Determinantenberechnung!)
r <- 3
# Anzahl der Sendeantennen variabel bis tmax
tmax <- 15
# Varianz der Kanalmatrix
sigma <- 1
# Varianz des Rauschens
N0 <- 1
# Varianz der Fehlermatrix
fsigma <- sigma/10
# Leistungsbeschränkung
P <- 10
# Anzahl der Durchläufe zur Durchschnittsermittlung
n <- 10000
# Initialisierung der Kapazitäten/Schranken
ev <- array(data = 0, dim = tmax)
sv <- array(data = 0, dim = tmax)
epu <- array(data = 0, dim = tmax)
epo <- array(data = 0, dim = tmax)
spu <- array(data = 0, dim = tmax)
spo <- array(data = 0, dim = tmax)
evk <- array(data = 0, dim = tmax)
svk <- array(data = 0, dim = tmax)
epuk <- array(data = 0, dim = tmax)
epok <- array(data = 0, dim = tmax)
spuk <- array(data = 0, dim = tmax)
spok <- array(data = 0, dim = tmax)
# Fehlertoleranz
fehler <- 0.000001

Hmax <- array(data = 0, dim = r)
Hi <- array(data = 0, dim = r)

# Einheitsmatrix
I <- diag(r)

# Schleife: Anzahl der Durchläufe für Mittelung
for (z in 1:n)
{
    # Schleife: variable Anzahl von Sendeantennen
    for (t in 1:tmax)
    {
        H <- matrix(nrow = r, ncol = t)
        for (i in 1:r)
        {
            for (j in 1:t)
            {
                H[i,j] <- complex(real = rnorm(1, mean=0, sd=sigma), imag = rnorm(1, mean=0, sd=sigma))
            }
        }
        HHstern <- H %*% Conj(t(H))

        # Senderseitig vollständig bekannt
        lambda <- eigen(HHstern)$values
        my <- P/min(r,t)
        for (i in 1:min(r,t))
        {        
            my <- my + lambda[i]^-1/min(r,t)
        }
        summe <- 0
        for (i in 1:min(r,t))
        {
            if (my>lambda[i]^-1) summe <- summe + (my - lambda[i]^-1)
        }
        if (abs(P - summe) > fehler)
        {
            my <- P/(min(r,t)-1)
            for (i in 1:(min(r,t)-1))
            {        
                my <- my + lambda[i]^-1/(min(r,t)-1)
            }
            summe <- 0
            for (i in 1:(min(r,t)-1))
            {
                if (my>lambda[i]^-1) summe <- summe + (my - lambda[i]^-1)
            }
            if (abs(P - summe) > fehler)
            {
                # Dieser Fall tritt so selten ein, dass er nicht berücksichtigt wird.
                cat("Fehler! Bitte starten Sie die Simulation neu. \n")
            }
        }
        summe <- 0
        for (i in 1:min(r,t))
        {
            if (my>0)
            if (log(my*lambda[i], base=exp(1))>0) summe <- summe + log(my*lambda[i], base=exp(1))
        }
        sv[t] <- sv[t] + summe

        # Empfängerseitig vollständig bekannt
        evmat <- I + P/t * HHstern
        evdet <- evmat[1,1]*evmat[2,2]*evmat[3,3] + evmat[2,1]*evmat[3,2]*evmat[1,3] + evmat[1,2]*evmat[2,3]*evmat[3,1] - evmat[1,3]*evmat[2,2]*evmat[3,1] - evmat[2,3]*evmat[3,2]*evmat[1,1] - evmat[1,2]*evmat[2,1]*evmat[3,3]
        # Imaginärteil = 0. Rundungsfehler müssen ausgeschlossen werden.
        evdet <- Re(evdet)
        ev[t] <- ev[t] + log(evdet, base=exp(1))

        # Empfängerseitig partiell bekannt (untere Schranke)
        epumat <- I + P/t/(P*sigma+N0) * HHstern
        epudet <- epumat[1,1]*epumat[2,2]*epumat[3,3] + epumat[2,1]*epumat[3,2]*epumat[1,3] + epumat[1,2]*epumat[2,3]*epumat[3,1] - epumat[1,3]*epumat[2,2]*epumat[3,1] - epumat[2,3]*epumat[3,2]*epumat[1,1] - epumat[1,2]*epumat[2,1]*epumat[3,3]
        # Imaginärteil = 0. Rundungsfehler müssen ausgeschlossen werden.
        epudet <- Re(epudet)
        epu[t] <- epu[t] + log(epudet, base=exp(1))

        # Empfängerseitig partiell bekannt (obere Schranke)
        epomat <- I + 1/sigma * HHstern
        epodet <- epomat[1,1]*epomat[2,2]*epomat[3,3] + epomat[2,1]*epomat[3,2]*epomat[1,3] + epomat[1,2]*epomat[2,3]*epomat[3,1] - epomat[1,3]*epomat[2,2]*epomat[3,1] - epomat[2,3]*epomat[3,2]*epomat[1,1] - epomat[1,2]*epomat[2,1]*epomat[3,3]
        # Imaginärteil = 0. Rundungsfehler müssen ausgeschlossen werden.
        epodet <- Re(epodet)
        epo[t] <- epo[t] + r*log(P*sigma/t/N0) + log(epodet, base=exp(1))

        # Senderseitig partiell bekannt (untere Schranke)
        spumat <- I +  P/t/(P*fsigma+N0) * HHstern
        spudet <- spumat[1,1]*spumat[2,2]*spumat[3,3] + spumat[2,1]*spumat[3,2]*spumat[1,3] + spumat[1,2]*spumat[2,3]*spumat[3,1] - spumat[1,3]*spumat[2,2]*spumat[3,1] - spumat[2,3]*spumat[3,2]*spumat[1,1] - spumat[1,2]*spumat[2,1]*spumat[3,3]
        # Imaginärteil = 0. Rundungsfehler müssen ausgeschlossen werden.
        spudet <- Re(spudet)
        spu[t] <- spu[t] + log(spudet, base=exp(1))

        # Senderseitig partiell bekannt (obere Schranke)
        for (i in 1:r)
        {
            Hmax[i] <- max(abs(H[i,1:t])^2)
            summeR <- 0
            if (t>1)
            for (l in 1:(t-1))
            {
                summeR <- summeR + abs(H[i,l])^2*sum(abs(H[i,(l+1):t])^2)
            }
            Hi[i] <- sqrt(summeR)
            spo[t] <- spo[t] + log(P*(4*Hi[i]+Hmax[i]+fsigma)/N0+1, base=exp(1))
        }

        # Hier folgen analoge Berechnungen für den Fall, dass P/t konstant gehalten wird - nicht P.
        # Dazu wird in den obigen Rechnungen P ersetzt durch t*P, so dass das fest definierte P die Leistungsbeschränkung pro Antenne angibt.

        # Senderseitig vollständig bekannt
        lambda <- eigen(HHstern)$values
        my <- P*t/min(r,t)
        for (i in 1:min(r,t))
        {        
            my <- my + lambda[i]^-1/min(r,t)
        }
        summe <- 0
        for (i in 1:min(r,t))
        {
            if (my>lambda[i]^-1) summe <- summe + (my - lambda[i]^-1)
        }
        if (abs(P*t - summe) > fehler)
        {
            my <- P*t/(min(r,t)-1)
            for (i in 1:(min(r,t)-1))
            {        
                my <- my + lambda[i]^-1/(min(r,t)-1)
            }
            summe <- 0
            for (i in 1:(min(r,t)-1))
            {
                if (my>lambda[i]^-1) summe <- summe + (my - lambda[i]^-1)
            }
            if (abs(P*t - summe) > fehler)
            {
                # Dieser Fall tritt so selten ein, dass er nicht berücksichtigt wird.
                print("Fehler! Bitte starten Sie die Simulation neu.")
            }
        }
        summe <- 0
        for (i in 1:min(r,t))
        {
            if (my>0)
            if (log(my*lambda[i], base=exp(1))>0) summe <- summe + log(my*lambda[i], base=exp(1))
        }
        svk[t] <- svk[t] + summe

        # Empfängerseitig vollständig bekannt
        evmat <- I + P * HHstern
        evdet <- evmat[1,1]*evmat[2,2]*evmat[3,3] + evmat[2,1]*evmat[3,2]*evmat[1,3] + evmat[1,2]*evmat[2,3]*evmat[3,1] - evmat[1,3]*evmat[2,2]*evmat[3,1] - evmat[2,3]*evmat[3,2]*evmat[1,1] - evmat[1,2]*evmat[2,1]*evmat[3,3]
        # Imaginärteil = 0. Rundungsfehler müssen ausgeschlossen werden.
        evdet <- Re(evdet)
        evk[t] <- evk[t] + log(evdet, base=exp(1))

        # Empfängerseitig partiell bekannt (untere Schranke)
        epumat <- I + P/(P*t*sigma+N0) * HHstern
        epudet <- epumat[1,1]*epumat[2,2]*epumat[3,3] + epumat[2,1]*epumat[3,2]*epumat[1,3] + epumat[1,2]*epumat[2,3]*epumat[3,1] - epumat[1,3]*epumat[2,2]*epumat[3,1] - epumat[2,3]*epumat[3,2]*epumat[1,1] - epumat[1,2]*epumat[2,1]*epumat[3,3]
        # Imaginärteil = 0. Rundungsfehler müssen ausgeschlossen werden.
        epudet <- Re(epudet)
        epuk[t] <- epuk[t] + log(epudet, base=exp(1))

        # Empfängerseitig partiell bekannt (obere Schranke)
        epomat <- I + 1/sigma * HHstern
        epodet <- epomat[1,1]*epomat[2,2]*epomat[3,3] + epomat[2,1]*epomat[3,2]*epomat[1,3] + epomat[1,2]*epomat[2,3]*epomat[3,1] - epomat[1,3]*epomat[2,2]*epomat[3,1] - epomat[2,3]*epomat[3,2]*epomat[1,1] - epomat[1,2]*epomat[2,1]*epomat[3,3]
        # Imaginärteil = 0. Rundungsfehler müssen ausgeschlossen werden.
        epodet <- Re(epodet)
        epok[t] <- epok[t] + r*log(P*sigma/N0) + log(epodet, base=exp(1))

        # Senderseitig partiell bekannt (untere Schranke)
        spumat <- I +  P/(P*t*fsigma+N0) * HHstern
        spudet <- spumat[1,1]*spumat[2,2]*spumat[3,3] + spumat[2,1]*spumat[3,2]*spumat[1,3] + spumat[1,2]*spumat[2,3]*spumat[3,1] - spumat[1,3]*spumat[2,2]*spumat[3,1] - spumat[2,3]*spumat[3,2]*spumat[1,1] - spumat[1,2]*spumat[2,1]*spumat[3,3]
        # Imaginärteil = 0. Rundungsfehler müssen ausgeschlossen werden.
        spudet <- Re(spudet)
        spuk[t] <- spuk[t] + log(spudet, base=exp(1))

        # Senderseitig partiell bekannt (obere Schranke)
        # Die Arrays Hmax und Hi wurden schon oben berechnet
        for (i in 1:r)
        {
            spok[t] <- spok[t] + log(P*t*(4*Hi[i]+Hmax[i]+fsigma)/N0+1, base=exp(1))
        }
    }
}

sv <- sv/n
ev <- ev/n
epu <- epu/n
epo <- epo/n
spu <- spu/n
spo <- spo/n
svk <- svk/n
evk <- evk/n
epuk <- epuk/n
epok <- epok/n
spuk <- spuk/n
spok <- spok/n

x11(width = 7, height = 7, pointsize = 12)
plot(1:15, 2*1:15, type="n", main="Verlauf der Kapazität bei wachsendem t und konst. P", xlab="t", ylab="C (Bits/Zeiteinheit)")
lines(sv, col="blue", lty="solid")
lines(ev, col="green", lty="dashed")
lines(epu, col="black", lty="dotted")
lines(epo, col="black", lty="dotted")
lines(spu, col="red", lty="dotdash")
lines(spo, col="red", lty="dotdash")

x11(width = 7, height = 7, pointsize = 12)
plot(1:15, 2*1:15, type="n", main="Verlauf der Kapazität bei wachsendem t und konst. P/t", xlab="t", ylab="C (Bits/Zeiteinheit)")
lines(svk, col="blue", lty="solid")
lines(evk, col="green", lty="dashed")
lines(epuk, col="black", lty="dotted")
lines(epok, col="black", lty="dotted")
lines(spuk, col="red", lty="dotdash")
lines(spok, col="red", lty="dotdash")

x11(width = 7, height = 7, pointsize = 12)
plot(1, 1, axes=FALSE, type="n", xlab="", ylab="", main="Legende")
lines(c(0,1,2), c(1.4,1.4,1.4), col="blue", lty="solid")
lines(c(0,1,2), c(1.2,1.2,1.2), col="green", lty="dashed")
lines(c(0,1,2), c(1,1,1), col="black", lty="dotted")
lines(c(0,1,2), c(0.8,0.8,0.8), col="red", lty="dotdash")
text(1, 1.3, "senderseitig vollständig bekannt")
text(1, 1.1, "empfängerseitig vollständig bekannt")
text(1, 0.9, "empfängerseitig partiell bekannt (obere und untere Schranke)")
text(1, 0.7, "senderseitig partiell bekannt (obere und untere Schranke)")