- S[j,mm,r] = -rho[mm,mm,r]*ps2[j,mm,r] + sum(phi[-j,mm,r] * Gram2[j,-j,r])
- if (abs(S[j,mm,r]) <= n*lambda*(pi[r]^gamma))
- phi[j,mm,r]=0
- else if(S[j,mm,r] > n*lambda*(pi[r]^gamma))
- phi[j,mm,r] = (n*lambda*(pi[r]^gamma)-S[j,mm,r]) / Gram2[j,j,r]
- else
- phi[j,mm,r] = -(n*lambda*(pi[r]^gamma)+S[j,mm,r]) / Gram2[j,j,r]
+ S[j, mm, r] <- -rho[mm, mm, r] * ps2[j, mm, r] +
+ sum(phi[-j, mm, r] * Gram2[j, -j, r])
+ if (abs(S[j, mm, r]) <= n * lambda * (pi[r]^gamma)) {
+ phi[j, mm, r] <- 0
+ } else if (S[j, mm, r] > n * lambda * (pi[r]^gamma)) {
+ phi[j, mm, r] <- (n * lambda * (pi[r]^gamma) - S[j, mm, r])/Gram2[j, j, r]
+ } else {
+ phi[j, mm, r] <- -(n * lambda * (pi[r]^gamma) + S[j, mm, r])/Gram2[j, j, r]
+ }