function  val = jsz_loadings_rho0cP(W, K1Q_X, rho0_cP, Sigma_cP, mats, dt)

  J = size(mats,2) ;
  N = size(K1Q_X,1) ;
  rho0d = 0 ;
  rho1d = ones(N,1) ;
  mats_periods = round(mats/dt) ;
  M = max(mats_periods) ;
  
  adjK1QX = jszAdjustK1QX(K1Q_X) ;
  K1Q_X = adjK1QX.K1Q_X ;
  m1 = adjK1QX.m1 ;
  
  K0Q_X = zeros(N,1) ;
  K0Q_X(m1) = 1  ;                 %% changed -- 1 instead of kinfQ
  loads_X_prelim = gaussian_loadings(mats_periods, K0Q_X, K1Q_X, zeros(N, N), rho0d*dt, rho1d*dt, dt) ;
  
  BX = loads_X_prelim.B ;
  alpha0_X = loads_X_prelim.A ;   %%% added
  WBXp = W * transpose(BX) ;           % N*N
  WBXpinv = inv(WBXp)     ;     %  N*N
  Sigma_X = WBXp\Sigma_cP /(WBXp')  ;  %% (W*BX')\Sigma.cP/(BX*W');
  
  % Now with Sigma_X in hand, compute loadings for AX
  loads_X = gaussian_loadings(mats_periods, K0Q_X, K1Q_X, Sigma_X, rho0d*dt, rho1d*dt, dt) ;
  
  AX1 = loads_X.A ;     %% 1*J         ### changed -- AX1 instead of AX
  
  alpha1_X = AX1 - alpha0_X;
  
  
%   ## Need to find what kinf should be to get the desired rho0_cP:
%   ## rt = 1'*Xt
%   ## cPt = (W*alpha0_X')*kinf + (W*alpha1_X') + (W*BX')*Xt
%   ## rt = 1'*(W*BX')^(-1)*[cPt -(W*alpha0_X')*kinf - (W*alpha1_X')]
%   ## --> rho0_cP = -1'*(W*BX')^(-1)*(W*alpha0_X')*kinf - 1'*(W*BX')^(-1)*(W*alpha1_X')
a0 = ones(1,N) / WBXp * W * alpha0_X'  ;
a1 = ones(1,N) / WBXp * W * alpha1_X'  ;
kinfQ = -(rho0_cP + a1)/a0;     %% -a0*kinf - a1 = rho0_cP
K0Q_X(m1) = kinfQ  ;

AX = alpha0_X*kinfQ + alpha1_X  ;

C = eye(J) - BX' * WBXpinv * W  ;
alpha0_cP = transpose(C * alpha0_X')  ;
alpha1_cP = transpose(C * alpha1_X')  ;

WAXp = W * AX' ;  %%  N*1
BcP = WBXpinv' * BX ;     %% N*J
AcP = AX * C'  ;

K1Q_cP = WBXp * K1Q_X * WBXpinv ;
K0Q_cP = WBXp * K0Q_X - K1Q_cP * WAXp ;

rho1_cP = WBXpinv' * ones(N,1) ;

val.AX=AX ; 
val.BX=BX ; 
val.AcP=AcP ;
val.BcP=BcP ;
val.K0Q_cP=K0Q_cP ; 
val.K1Q_cP=K1Q_cP ;
val.rho0_cP=rho0_cP ;
val.K0Q_X=K0Q_X ;
val.K1Q_X=K1Q_X ;
val.rho1_cP=rho1_cP ;
val.alpha0_X=alpha0_X ;
val.alpha1_X=alpha1_X ;
val.alpha0_cP=alpha0_cP ;
val.alpha1_cP=alpha1_cP ;
val.m1=m1 ;

end


