function [ W_ind, W_dat ] = medfuncSystemMatrix( theta, N, P_num, delta )
%UNTITLED Summary of this function goes here
% Detailed explanation goes here
N2 = N^2;
M = length(theta)*P_num;
W_ind = zeros(M, 2*N);
W_dat = zeros(M, 2*N);
t = (-(P_num-1)/2:(P_num-1)/2)*delta;
if N<=10 && length(theta)<=5
x = (-N/2:1:N/2)*delta;
y = (-N/2:1:N/2)*delta;
plot(x,meshgrid(y,x),'k');
hold on;
plot(x,meshgrid(x,y),y,'k');
axis([-N/2-5, N/2+5, -N/2-5, N/2+5]);
text(0,-0.4*delta, '0');
end
for jj=1:length(theta)
for ii=1:1:P_num
u = zeros(1,2*N);
v = zeros(1,2*N);
th = theta(jj);
if th >= 180 || th<0
error('输入角度有问题');
elseif th==90
if N<=10 && length(theta)<=5
xx = (-N/2-2:0.01:N/2+2)*delta;
yy = t(ii);
plot(xx,yy,'b');
hold on
end
if t(ii)>=N/2*delta || t(ii)<=-N/2*delta
continue;
end
kout = N*ceil(N/2 - t(ii)/delta);
kk = (kout - (N-1)):1:kout;
u(1:N) = kk;
v(1:N) = ones(1,N)*delta;
elseif th==0
if N<=10 && length(theta)<=5
yy = (-N/2-2:0.01:N/2+2)*delta;
xx = t(ii);
plot(xx,yy,'b');
hold on
end
if t(ii)>=N/2*delta || t(ii)<=-N/2*delta
continue;
end
kin = ceil(N/2+t(ii)/delta);
kk = kin:N:(kin+N*(N-1));
u(1:N) = kk;
v(1:N) = ones(1,N)*delta;
else
if th>90
th_temp = th-90;
elseif th<90
th_temp = 90-th;
end
th_temp = th_temp*pi/180;
b = t/cos(th_temp);
m = tan(th_temp);
y1d = -N/2*delta*m + b(ii);
y2d = N/2*delta*m + b(ii);
if N<=10 && length(theta)<=5
xx = (-N/2-2:0.01:N/2+2)*delta;
if th<90
yy = -m*xx + b(ii);
elseif th>90
yy = m*xx + b(ii);
end
plot(xx,yy,'b');
hold on;
end
if y1d<-N/2*delta & y2d<-N/2*delta
continue;
end
if y1d>N/2*delta & y2d>-N/2*delta
continue;
end
if y1d<=N/2*delta & y1d>=-N/2*delta & y2d>N/2*delta
yin = y1d;
d1 = yin - floor(yin/delta)*delta;
kin = N*floor(N/2 - yin/delta) + 1;
yout = N/2*delta;
xout = (yout-b(ii))/m;
kout = ceil(xout/delta) + N/2;
elseif y1d<=N/2*delta & y1d >= -N/2*delta & y2d >= -N/2*delta & y2d<N/2*delta
yin = y1d;
d1 = yin - floor(yin/delta)*delta;
kin = N*floor(N/2-yin/delta) + 1;
yout = y2d;
kout = N * floor(N/2 - yout/delta) + N;
elseif y1d<-N/2*delta & y2d>N/2*delta
yin = -N/2*delta;
xin = (yin-b(ii))/m;
d1 = N/2*delta + (floor(xin/delta)*delta*m + b(ii));
kin = N*(N-1) + N/2 + ceil(xin/delta);
yout = N/2*delta;
xout = (yout - b(ii))/m;
kout = ceil(xout/delta) + N/2;
elseif y1d<-N/2*delta & y2d>=-N/2*delta & y2d<N/2*delta
yin = -N/2*delta;
xin = (yin-b(ii))/m;
d1 = N/2*delta + (floor(xin/delta)*delta*m + b(ii));
kin = N*(N-1) + N/2 + ceil(xin/delta);
yout = y2d;
kout = N*floor(N/2 - yout/delta) + N;
else
continue;
end
k = kin;
c=0;
d2 = d1 + m*delta;
while k>=1 && k<=N2
c = c + 1;
if d1>=0 && d2 >delta
u(c) = k;
v(c) = (delta-d1)*sqrt(m^2+1)/m;
if k>N && k~=kout
k = k - N;
d1 = d1 - delta;
d2 = d1 + m*delta;
else
break;
end
elseif d1>=0 && d2 == delta
u(c) = k;
v(c) = delta*sqrt(m^2 + 1);
if k>N && k~=kout
k = k-N+1;
d1 = 0;
d2 = d1 + m*delta;
else
break;
end
elseif d1>=0 && d2<delta
u(c) = k;
v(c) = delta*sqrt(m^2+1);
if k~=kout
k = k+1;
d1 = d2;
d2 = d1 + m*delta;
else
break;
end
elseif d1<=0 && d2 >=0 && d2<=delta
u(c) = k;
v(c) = d2*sqrt(m^2+1)/m;
if k~=kout
k = k+1;
d1 = d2;
d2 = d1 + m*delta;
else
break;
end
elseif d1<=0 && d2>delta
u(c) = k;
v(c) = delta*sqrt(m^2+1)/m;
if k>N && k~=kout
k = k - N;
d1 = -delta + d1;
d2 = d1 + m*delta;
else
break;
end
end
end
if th < 90
u_temp = zeros(1,2*N);
if any(u) == 0
continue;
end
ind = u>0;
for k=1:length(u(ind))
r = rem(u(k),N);
if r==0
u_temp(k) = u(k) - N + 1;
else
u_temp(k) = u(k) - 2*r + N + 1;
end
end
u = u_temp;
end
end
W_ind((jj-1)*P_num+ii,:) = u;
W_dat((jj-1)*P_num+ii,:) = v;
end
end
