% Daten :  n = Phantomgroesse
%               Q = Quellpunkt, gegeben in Polarkoordinaten (r,phi)
%               alpha=Auslenkungen vom Zentralstrahl
 
% n=10;
% phi=0.678;
% r=2*n;


% Ausgangsintensitaet% Matrix einlesen?

 rmat = input('Matrix einlesen (j/n)? ', 's');
 if rmat == 'j'
    fn = input('Dateiname : ','s');
    load(fn);
    SC = M;
    n = size(SC,1)/2; 
    SM = reshape(SC',4*n*n,1);
 else
    n=input('Gittergroesse (4xn^2) ');
    SC = reshape(make_scross(n),2*n,2*n);
    SM = reshape(SC',4*n*n,1);
 end

% Eingabe der spez. Werte

 phi=input('Quellpunkt in Polarkoordinaten (phi) ');
 r=2*n;
% alpha=input('Streubreite (Auslenkung) ')

 alpha=[-4:4]*pi/48;
 Q=r*[cos(phi);sin(phi)];
 
 in = input('Anfangsintensitaet (def Iin=100) :');
 if size(in, 1) == 0, Iin  = 100, else Iin = in, end;

% Schrittweite step
 step = input('Schrittweite :');
 

% Berechnen der Schnittpunkte

I = [];
L = sparse(step*length(alpha), 4*n*n);
eq = 0;

b_break=0;

d2 = sqrt(2);
beta2 = pi/4;


FIG = figure('Name', 'CT Simulation');
%hold off,
set(gca, 'YDir', 'normal');
set(gca, 'DrawMode', 'fast');
tic;
for k=1:step
   phi=k*2*pi/step;
   Q = r*[cos(phi);sin(phi)];
   d = sqrt((Q(1)-1)^2+(Q(2)-1)^2);
   d_alpha = 0% asin(sin(pi/2-phi+beta2)*d2/d);
   s_a=tan(phi+alpha+d_alpha); 

   % Zeichne Phantom; besteht aus regelmaessigem Gitter
   % mit Maschenweite 1 von -n bis +n

   [x,y]=meshgrid([-n:1:n],[-n:1:n]);
 
   % SC = reshape(make_scross(n),2*n,2*n);
%  hold off,
   colormap(gray);
   image(0.5-n,0.5-n,(SC'+0.5)*60); 
   hold on, 

   plot(x,y,'-g'), 
   hold on, 
   plot(y,x,'-g'),

   % Zeichne Kreis, auf dem sich die Quelle bewegt, und Quelle

   kreis_x=r*cos(2*pi*[0:200]/200);
   kreis_y=r*sin(2*pi*[0:200]/200);
   plot(kreis_x,kreis_y,'-r')                           
   plot(Q(1),Q(2),'or')  

   % Zeichne Zentralstrahl
   %plot([Q(1) -Q(1)],[Q(2) -Q(2)],'-y') 
   % s_a=tan(phi+alpha);
   P=Q*ones(1,length(alpha))-[cos(alpha+phi+d_alpha);sin(alpha+phi+d_alpha)] ...
   *diag(cos(alpha))*2*r;
   for j=1:length(alpha),
       plot([Q(1) P(1,j)],[Q(2) P(2,j)],'-c')
   end
   plot([Q(1) P(1,5)],[Q(2) P(2,5)],'-r')
   axis([-r r -r r])
   axis('square')
   pause(1),
   clf,

b_break=0;
   
   spunkt

   % Austrittsintensitaeten

   ST = reshape(SM(B), 4*n+1, 9);
   Iout = exp(sum(-ST.*len))*Iin;
   I=[I; (log(Iout/Iin))'];

   [row, col] = find(len);

   for c=1:length(alpha)
     eq=eq+1;
     va = row(col==c);   
     inx = B(va, c);
     L(eq,inx)=-len(va,c)';
     if norm(len(va,c))==0, b_break=1; break , end
     if abs(log(Iout(c)/Iin)-L(eq,:)*SM)>10^(-6), b_break=1; break , end
   end
   if b_break==1, break, end
end
toc
figure('Name', 'Koeffizientenverteilung')
spy(L);



figure('Name', 'Berechnetes Phantom')
colormap(gray);

tic;
% TT = pinv(full(L),10^(-12))*I;
TT = L\I;
toc
T = reshape(TT, r, r);
image(0.5-n,0.5-n,(T+0.5)*60); 
hold on, set(gca, 'YDir', 'normal');
axis([-n n -n n])
axis('square')

I_g=I+randn(size(I))*6*10^(-2);
TT_g = L\I_g;
T_g = reshape(TT_g, r, r);
figure
image(0.5-n,0.5-n,(T_g+0.5)*60); 
hold on, set(gca, 'YDir', 'normal');
axis([-n n -n n])
axis('square')
