close all
clc
d0=23.7;
ci=0.075;
cc=1.16*ci;
cq=cc;


r=d0/2;                                 %Radius der Austrittsdüse
phi=(0:pi/8:2*pi);                     %Schwerpunktlinien
alpha=(0+pi/16:pi/8:2*pi);             %Raster 


%Schwerpunktsberechnung Elementardüsen
r1=12.5;
r2=8.5;
r3=5.5;
r4=3.5;
r5=2.5;
[x_Kreis1, y_Kreis1]=pol2cart(alpha,r1);
[x_Kreis2, y_Kreis2]=pol2cart(alpha,r2);
[x_Kreis3, y_Kreis3]=pol2cart(alpha,r3);
[x_Kreis4, y_Kreis4]=pol2cart(alpha,r4);
[x_Kreis5, y_Kreis5]=pol2cart(alpha,r5);

%Längere Seite des Trapez ---> Seitea
for i=1:length(alpha)
    for j=i:length(alpha)
        Seitea= sqrt((x_Kreis1(i)-x_Kreis1(j))^2+(y_Kreis1(i)-y_Kreis1(j))^2);
    end 
end
%Kürzere Seite des Trapez --->Seiteb
for i=1:length(alpha)
    for j=i:length(alpha)
        Seiteb=sqrt((x_Kreis2(i)-x_Kreis2(j))^2+(y_Kreis2(i)-y_Kreis2(j))^2);
    end
end






figure('Name','Düse','NumberTitle','off')
hold on
plot(x_Kreisa,y_Kreis1,'r')
plot(MR(:,1),MR(:,2),'+','Markersize',3)
grid on

%Stromdichteverteilung

%Brennergeometrie liegt in x,y-Ebene

Z=100;                %Laufweite 


 X=-300:10:300;
 Y=-300:10:300;



a=Raster_R;
b=a;
xj=MR(:,1);
yj=MR(:,2);

for i=1:length(X)
    for j=1:length(Y)
        ImpulsD(i,j)=sum(1/((4*ci*ci*(X(i)-xj).^2)/(Raster_R.^2))*exp(-((Y(j)-yj).^2)/(ci*ci*(X(i)-xj).^2))) ;
    end
end



figure('Name','Plot')
surf(X,Y,ImpulsD)