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

phi = linspace(0,2*pi);                 %Winkel zur Kreisdarstellung
R=275;                                  %Radius Düsenanordnung
r=d0/2;                                 %Radius der Austrittsdüse

theta1=0;                               %Düsenposition (Winkel) im polaren KS auf R
theta2=theta1+2*pi/8;
theta3=theta2+2*pi/8;
theta4=theta3+2*pi/8;
theta5=theta4+2*pi/8;
theta6=theta5+2*pi/8;
theta7=theta6+2*pi/8;
theta8=theta7+2*pi/8;

di=200;                                 %Innendurchmesser Luft
da=416;                                 %Außendurchmesser Luft
ri=di/2;                                    
ra=da/2;

%---------------------------------------
%Sekundärgasdüsen positionieren
%---------------------------------------
[X,Y]=pol2cart(phi,R);
[x1,y1]=pol2cart(theta1,R);
[x2,y2]=pol2cart(theta2,R);
[x3,y3]=pol2cart(theta3,R);
[x4,y4]=pol2cart(theta4,R);
[x5,y5]=pol2cart(theta5,R);
[x6,y6]=pol2cart(theta6,R);
[x7,y7]=pol2cart(theta7,R);
[x8,y8]=pol2cart(theta8,R);
M_Sek=[x1,y1;x2,y2;x3,y3;x4,y4;x5,y5;x6,y6;x7,y7;x8,y8];
x_Sek=M_Sek(:,1);
y_sek=M_Sek(:,2);

% Bestimmen der Punkte eines bestimmten Rasters, die innerhalb des Kreises liegen 
Raster = 0.2; 
x = -(r+Raster):Raster:r+Raster; 
y = -(r+Raster):Raster:r+Raster;  
n=1; 
for i=1:length(x) 
    for j=1:length(y)
        for k=1:length(x_Sek)
            if x(i)^2+y(j)^2 <= r^2 
                M_Ed(n,1:2) = [x(i)+M_Sek(k,1),y(j)+M_Sek(k,2)]; 
                n = n+1;
            end
        end 
    end 
end 


% Bestimmen der Punkte eines bestimmten Rasters, die innerhalb des Ringspaltes liegen 
ri=100;
ra=205;
[x_Kreisi y_Kreisi]=pol2cart(phi,ri);
[x_Kreisa y_Kreisa]=pol2cart(phi,ra);
Raster_R=5;
xR = -(ra+Raster_R):Raster_R:ra+Raster_R; 
yR = -(ra+Raster_R):Raster_R:ra+Raster_R; 
n=1; 
for i=1:length(xR) 
    for j=1:length(yR) 
        if xR(i)^2+yR(j)^2 <= ra^2
            if xR(i)^2+yR(j)^2 >= ri^2
                MR(n,1:2) = [xR(i),yR(j)]; 
                n = n+1;
            end
        end 
    end 
end

figure('Name','Brennergeometrie','NumberTitle','off')
hold on
plot(M_Ed(:,1),M_Ed(:,2),'+','Markersize',3)
plot(x_Kreisi,y_Kreisi,'r')
plot(x_Kreisa,y_Kreisa,'r')
plot(MR(:,1),MR(:,2),'+','Markersize',3)
grid on
%----------------------------------
%Stromdichteverteilung
%----------------------------------
%Brennergeometrie liegt in x,y-Ebene

Z=3000;                %Laufweite 


[X,Y]=meshgrid(-300:1:300,-300:1:300);
a=Raster;
b=a;
xmj=M_Ed(:,1);
ymj=M_Ed(:,2);
zmj=0;
SDViSek=sum(0.25*(erf(((Y-ymj)+a/2)./(ci*(Z-zmj)))-erf(((Y-ymj)-a/2)./(ci*(Z-zmj)))).*(erf(((X-xmj)+b/2)./(ci*(Z-zmj)))-erf(((X-xmj)-b/2)./(ci*(Z-zmj)))));


%Ringspalt

a=Raster_R;
b=a;
xmjR=MR(:,1);
ymjR=MR(:,2);
zmjR=0;
SDViLuft=sum(0.25*(erf(((Y-ymjR)+a/2)./(ci*(Z-zmjR)))-erf(((Y-ymjR)-a/2)./(ci*(Z-zmjR)))).*(erf(((X-xmjR)+b/2)./(ci*(Z-zmjR)))-erf(((X-xmjR)-b/2)./(ci*(Z-zmjR)))))


mesh(SDViSek,X,Y)