Browse Source

Gitlab Upgrade

clayshieh 3 years ago
commit
bf5c09911f
21 changed files with 3782 additions and 0 deletions
  1. 44
    0
      Inputs.csv
  2. 217
    0
      Modified MCode/SN_per_bench_solver.m
  3. 163
    0
      Modified MCode/periodic_benchmark.m
  4. 382
    0
      Vanilla MCode/New_SN_bench_solver_MC.m
  5. 212
    0
      Vanilla MCode/SN_hom.m
  6. 207
    0
      Vanilla MCode/SN_per_bench_solver.m
  7. 143
    0
      Vanilla MCode/periodic_benchmark.m
  8. 32
    0
      compare.py
  9. 202
    0
      csv/A.csv
  10. 202
    0
      csv/B.csv
  11. 102
    0
      csv/L.csv
  12. 101
    0
      csv/SCAL.csv
  13. 101
    0
      csv/SF.csv
  14. 200
    0
      csv/X.csv
  15. 202
    0
      csv/Y.csv
  16. 202
    0
      csv/Z.csv
  17. 100
    0
      csv/extra.csv
  18. 101
    0
      csv/h.csv
  19. 21
    0
      csv/xsn.csv
  20. 10
    0
      steady-state.config
  21. 838
    0
      steady-state.py

+ 44
- 0
Inputs.csv
File diff suppressed because it is too large
View File


+ 217
- 0
Modified MCode/SN_per_bench_solver.m View File

@@ -0,0 +1,217 @@
1
+%THIS IS THE AUXILIARY FUNCTION OF THE MAIN CODE "periodic_benchmark"
2
+
3
+%RICHARD VASQUES & NITIN KUMAR YADAV
4
+
5
+
6
+function[Z,extra,n1,B,L,A,h] = SN_per_bench_solver(T,m1,m2,n,N,Es1,Es2,Et1,Et2,yo,y_,Q1,Q2,u,wt,a)
7
+
8
+%Periodic Medium........
9
+interval=T/n;
10
+m12=zeros(2,1);
11
+m12(1)=m1;m12(2)=m2;
12
+mm=0;
13
+i=0;
14
+s=0;
15
+if a>1
16
+    if a<m2/interval+2
17
+		i=i+1;
18
+		x1=(a-1)*interval;
19
+		s=s+x1;
20
+		x(i,1)=s;
21
+		x(i,2)=2;
22
+	else
23
+    	i=i+1;
24
+		x1=(a-m2/interval-1)*interval;
25
+		s=s+x1;
26
+		x(i,1)=s;
27
+		x(i,2)=mod(mm,2)+1;
28
+       mm=mm+1;
29
+    end
30
+end
31
+
32
+while s<T
33
+    i=i+1;
34
+    x1=m12(mod(mm,2)+1);
35
+    s=s+x1;
36
+    if s<=T
37
+        x(i,1)=s;
38
+    else
39
+        x(i,1)=T;
40
+    end
41
+    x(i,2)=mod(mm,2)+1;
42
+    mm=mm+1;
43
+end
44
+
45
+H=T/n;n1=1;i=1;j=1;
46
+extra=zeros(n,1);
47
+t1=i*H; %%%%%%WHATS THE POINT OF THIS???
48
+dlmwrite('xsn_m.txt',x,'delimiter',',','precision',18);
49
+display('saved xsn');
50
+L(1,1)=0;L(1,2)=x(j,2);
51
+if t1==x(j,1)
52
+    display('hi');
53
+    extra(i)=1;
54
+    h(n1)=x(j,1)/2;
55
+    L(n1+1,1)=h(n1);
56
+    L(n1+1,2)=x(j,2);
57
+    h(n1+1)=h(n1);
58
+    L(n1+2,1)=x(j,1);
59
+    L(n1+2,2)=3;
60
+    n1=n1+2;
61
+     i=i+1;
62
+     t1=i*H;
63
+else
64
+    while t1<=x(j,1)
65
+        h(n1)=H;
66
+        display(H);
67
+        L(n1+1,1)=i*H;
68
+        if L(n1+1,1)==x(j,1)                 % checking is point is interface point?
69
+            L(n1+1,2)=3;
70
+        else
71
+            L(n1+1,2)=x(j,2);
72
+        end
73
+        n1=n1+1;
74
+        i=i+1;
75
+        t1=i*H;
76
+    end
77
+end
78
+j=2;
79
+if x(1,1)~=T
80
+    while i<=n
81
+        while t1<=x(j,1)
82
+                L(n1+1,1)=i*H;
83
+                if L(n1+1,1)==x(j,1)          %checking is point is interface point
84
+                    L(n1+1,2)=3;
85
+                else
86
+                    L(n1+1,2)=x(j,2);
87
+                end
88
+                h(n1)=L(n1+1,1)-L(n1,1);
89
+                n1=n1+1;
90
+                i=i+1;
91
+                t1=i*H;
92
+        end
93
+        j=j+1;
94
+        if L(n1,1)==T && L(n1,2)==3 && L(n1-1,2)==3
95
+            display('boop');
96
+            display(n1);
97
+            i=i-1; 
98
+            extra(i)=extra(i)+1;
99
+            L(n1,1)=(x(j-1,1)+x(j-2,1))/2;
100
+            L(n1,2)=x(j-1,2);
101
+            h(n1-1)=L(n1,1)-L(n1-1,1);
102
+            n1=n1+1;
103
+            display(n1);
104
+            L(n1,1)=x(j-1,1);
105
+            L(n1,2)=3;
106
+            
107
+            h(n1-1)=L(n1,1)-L(n1-1,1);
108
+            i=i+1;
109
+        end       
110
+    end
111
+end
112
+n1=n1-1;
113
+A=zeros(N*n1,N*n1); B=zeros(N*n1,1);
114
+Et12=zeros(2,1); Es12=zeros(2,1); Q12=zeros(2,1);
115
+Et12(1)=Et1;Et12(2)=Et2;
116
+Es12(1)=Es1;Es12(2)=Es2;
117
+Q12(1)=Q1;Q12(2)=Q2;
118
+
119
+
120
+% Diagonal Block of matrix up to N/2.................................
121
+for t=1:N/2
122
+    s=(t-1)*n1;
123
+    A(s+1,s+1)=-u(t)*(1/h(1))+Et12(L(1,2))-Es12(L(1,2))*wt(t)/2;
124
+    A(s+1,s+2)=u(t)*(1/h(1));
125
+    B(s+1)=Q12(L(1,2))/2;
126
+    for i=2:n1-1
127
+        if L(i,2)==3
128
+            if i==n1-1
129
+                A(s+i,s+i)=-u(t)*(1/h(i)+1/(h(i)+h(i+1)))+Et12(L(i+1,2))-Es12(L(i+1,2))*wt(t)/2;
130
+                A(s+i,s+i+1)=u(t)*(1/h(i)+1/h(i+1));
131
+                B(s+i)=u(t)*y_*(h(i)/(h(i+1)*(h(i)+h(i+1))))+Q12(L(i+1,2))/2;
132
+            else
133
+                A(s+i,s+i)=-u(t)*(1/h(i)+1/(h(i)+h(i+1)))+Et12(L(i+1,2))-Es12(L(i+1,2))*wt(t)/2;
134
+                A(s+i,s+i+1)=u(t)*(1/h(i)+1/h(i+1));
135
+                A(s+i,s+i+2)=-u(t)*(h(i)/(h(i+1)*(h(i)+h(i+1))));
136
+                B(s+i)=Q12(L(i+1,2))/2;
137
+        
138
+            end
139
+        else
140
+            A(s+i,s+i-1)=-u(t)*h(i)/(h(i-1)*(h(i-1)+h(i)));
141
+            A(s+i,s+i)=u(t)*(h(i)-h(i-1))/(h(i)*h(i-1))+Et12(L(i,2))-Es12(L(i,2))*wt(t)/2;
142
+            A(s+i,s+i+1)=u(t)*h(i-1)/(h(i)*(h(i-1)+h(i)));
143
+            	B(s+i)=Q12(L(i,2))/2;
144
+        end
145
+    end
146
+    A(s+n1,s+n1-1)=-u(t)*h(n1)/(h(n1-1)*(h(n1-1)+h(n1)));
147
+    A(s+n1,s+n1)=u(t)*(h(n1)-h(n1-1))/(h(n1)*h(n1-1))+Et12(L(n1,2))-Es12(L(n1,2))*wt(t)/2;
148
+    B(s+n1)=-u(t)*y_*h(n1-1)/(h(n1)*(h(n1-1)+h(n1)))+Q12(L(n1,2))/2;
149
+    
150
+    % Blocks from N/2 to N........................................
151
+    a=0;
152
+    for p=1:N/2
153
+        S=(N/2+p-1)*n1;
154
+        for i=2:n1
155
+            if L(i,2)==3
156
+                A(s+i,S+i-1)=-Es12(L(i+1,2))*wt(N/2+p)/2;
157
+            else
158
+                A(s+i,S+i-1)=-Es12(L(i,2))*wt(N/2+p)/2;
159
+            end
160
+        end
161
+        a=a+(Es12(L(1,2))*wt(N/2+p)*yo/2);
162
+    end
163
+    B(s+1)=B(s+1)+a;
164
+end
165
+    
166
+% Diagonal Block of matrix from N/2+1 to N.........................
167
+for t=N/2+1:N
168
+    s=(t-1)*n1;
169
+    A(s+1,s+1)=u(t)*(h(2)-h(1))/(h(2)*h(1))+Et12(L(1,2))-Es12(L(1,2))*wt(t)/2;
170
+    A(s+1,s+1+1) = u(t)*h(1)/(h(2)*(h(1)+h(2)));
171
+    B(s+1)=u(t)*yo*h(2)/(h(1)*(h(1)+h(2)))+Q12(L(1,2))/2;
172
+    for i=2:n1-1
173
+        if L(i+1,2)==3
174
+            if i==2
175
+                A(s+i,s+i)=u(t)*(1/h(i)+1/(h(i)+h(i-1)))+Et12(L(i,2))-Es12(L(i,2))*wt(t)/2;
176
+                A(s+i,s+i-1)=-u(t)*(1/h(i)+1/h(i-1));
177
+                B(s+i)=-u(t)*yo*(h(i)/(h(i-1)*(h(i)+h(i-1))))+Q12(L(i,2))/2;
178
+            else
179
+                A(s+i,s+i)=u(t)*(1/h(i)+1/(h(i)+h(i-1)))+Et12(L(i,2))-Es12(L(i,2))*wt(t)/2;
180
+                A(s+i,s+i-1)=-u(t)*(1/h(i)+1/h(i-1));
181
+                A(s+i,s+i-2)=u(t)*(h(i)/(h(i-1)*(h(i)+h(i-1))));
182
+                B(s+i)=Q12(L(i,2))/2;
183
+                if s+i ==105
184
+                    fprintf('105');
185
+                end
186
+	        end
187
+        else
188
+            A(s+i,s+i-1)=-u(t)*h(i+1)/(h(i)*(h(i)+h(i+1)));
189
+            A(s+i,s+i)=u(t)*(h(i+1)-h(i))/(h(i+1)*h(i))+Et12(L(i+1,2))-Es12(L(i+1,2))*wt(t)/2;
190
+            A(s+i,s+i+1)=u(t)*h(i)/(h(i+1)*(h(i)+h(i+1)));
191
+            B(s+i)=Q12(L(i+1,2))/2;
192
+	    end
193
+    end
194
+    A(s+n1,s+n1)=u(t)*(1/h(n1))+Et12(L(n1,2))-Es12(L(n1,2))*wt(t)/2;
195
+    A(s+n1,s+n1-1)=-u(t)*(1/h(n1));
196
+    B(s+n1)=Q12(L(n1,2))/2;
197
+    
198
+    % Blocks from 1 to N/2........................................
199
+    a=0;
200
+    for p=1:N/2
201
+        S=(p-1)*n1;
202
+        for i=1:n1-1
203
+            if L(i+1,2)==3
204
+                A(s+i,S+i+1)=-Es12(L(i,2))*wt(p)/2;
205
+            else
206
+                A(s+i,S+i+1)=-Es12(L(i+1,2))*wt(p)/2;
207
+            end
208
+        end
209
+        a=a+(Es12(L(n1,2))*wt(p)*y_/2);
210
+    end
211
+    B(s+n1)=B(s+n1)+a;
212
+end
213
+
214
+Z=A\B; % Solving for angular flux.
215
+
216
+return
217
+end

+ 163
- 0
Modified MCode/periodic_benchmark.m View File

@@ -0,0 +1,163 @@
1
+%THIS CODE GENERATES ENSEMBLE-AVERAGED BENCHMARK RESULTS FOR THE
2
+%STEADY-STATE, MONOENERGETIC TRANSPORT EQUATION IN A PERIODIC MEDIUM
3
+%COMPOSED OF TWO MATERIALS WITH ISOTROPIC SCATTERING AND ISOTROPIC SOURCE
4
+%IN ROD GEOMETRY. THE TOTAL NUMBER OF REALIZATIONS IS DEFINED BY THE
5
+%SPATIAL DISCRETIZATION.
6
+
7
+%THIS CODE CALLS THE SOLVER ROUTINE "SN_per_bench_solver", USING 2-POINT CENTRAL
8
+%DIFFERENCING (ORDER H^2).
9
+
10
+clear all
11
+clc
12
+
13
+%------------------------------------------------------------
14
+%INPUTS
15
+T= 20;       % total length of the system
16
+n= 100;     % # of points (discretization)
17
+
18
+yo=0.;      % left boundary condition
19
+y_=0.;      % right boundary condition
20
+
21
+m1= 1.0;    % thickness of material 1 layers
22
+m2= 1.0;    % thickness of material 2 layers
23
+
24
+Et1= 1.0;    % total cross section Sigma_t1 of material 1
25
+cc= 0.5;    % Scattering Ratio c1 of material 1
26
+Q1= 1.0;     % homogeneous isotropic Source of material 1
27
+
28
+Et2=0.;      % total cross section Sigma_t2 of material 2
29
+cc2= 0.;    % Scattering Ratio c2 of material 2
30
+Q2= 0.;     % homogeneous isotropic Source of material 2
31
+%------------------------------------------------------------
32
+%------------------------------------------------------------
33
+
34
+Es1=cc*Et1;
35
+Ea1= Et1-Es1;
36
+Es2=cc2*Et2;
37
+Ea2= Et2-Es2;
38
+
39
+%weights and directions
40
+N=2;
41
+wt(1)= 1;
42
+wt(2)=1;
43
+u(1)=-1;
44
+u(2)=1;
45
+
46
+a=1;
47
+reflec=0;reflec2=0;
48
+transm=0;transm2=0;
49
+SF=zeros(n+1,1);
50
+SF2=zeros(n+1,1);
51
+cond=0;
52
+
53
+%MAIN LOOP
54
+total = (m1+m2)/(T/n);
55
+
56
+while(a<(total+1))
57
+    
58
+    fprintf('Problem %d of %d\n', a, total);
59
+    [Z,extra,n1,B,L,A,h] = SN_per_bench_solver(T,m1,m2,n,N,Es1,Es2,Et1,Et2,yo,y_,Q1,Q2,u,wt,a);
60
+    dlmwrite('L_m.txt',L,'delimiter',',','precision',18);
61
+    dlmwrite('Z_m.txt',Z,'delimiter',',','precision',18);
62
+    dlmwrite('A_m.txt',A,'delimiter',',','precision',18);
63
+    dlmwrite('B_m.txt',B,'delimiter',',','precision',18);
64
+    dlmwrite('extra_m.txt',extra,'delimiter',',','precision',18);
65
+    dlmwrite('h_m.txt',transpose(h),'delimiter',',','precision',18);
66
+    
67
+            % Adjusting points..........................
68
+            X=zeros(n*N,1);
69
+            for i=1:N/2
70
+                X((i-1)*n+1)=Z((i-1)*n1+1);
71
+                k=2;
72
+                for j=2:n
73
+                    %display(k);
74
+                    k=k+extra(j-1);
75
+                    %display(k);
76
+                    X((i-1)*n+j)=Z((i-1)*n1+k);
77
+                    k=k+1;
78
+                end
79
+            end            
80
+            for i=N/2+1:N
81
+                k=1;
82
+                for j=1:n
83
+                    %display(k);
84
+                    k=k+extra(j);
85
+                    %display((i-1)*n1+k);
86
+                    X((i-1)*n+j)=Z((i-1)*n1+k);
87
+%                     if (i-1)*n+j == 102
88
+%                         display('here');
89
+%                         display((i-1)*n1+k);
90
+%                         display(Z((i-1)*n1+k));
91
+%                     end
92
+                    k=k+1;
93
+                end
94
+            end
95
+            
96
+            Y=zeros(N*(n+1),1);
97
+            % Adding boundary conditions...............
98
+            i=1;
99
+            for t=1:N/2
100
+                s=(t-1);
101
+                for j=i:t*n+s
102
+                    Y(j)=X(j-s);
103
+                end
104
+                Y(j+1)=y_;
105
+                i=j+2;
106
+            end
107
+            i=i-1;
108
+            for t=N/2+1:N
109
+                Y(i+1)=yo;
110
+                i=i+1;
111
+                for j=i+1:t*n+t
112
+                    Y(j)=X(j-t);
113
+                end 
114
+                i=j;
115
+            end
116
+            
117
+            dlmwrite('X_m.txt',X,'delimiter',',','precision',18);
118
+            dlmwrite('Y_m.txt',Y,'delimiter',',','precision',18);
119
+            
120
+            % Calculating Reflection, Transmission and Scalar Flux......     
121
+            RL=0;
122
+            for t=1:N/2
123
+                s=(t-1)*n;
124
+                RL=RL+abs(wt(t)*u(t)*Y(s+t));
125
+            end
126
+            TR=0;
127
+            for t=N/2+1:N
128
+                s=(t-1)*n;
129
+                TR=TR+wt(t)*u(t)*Y(s+n+t);
130
+            end
131
+            m=n+1;
132
+            SCAL=zeros(m,1);
133
+            for t=1:N
134
+                s=(t-1)*m;
135
+                for i=1:m
136
+                    SCAL(i)=SCAL(i)+wt(t)*Y(s+i);
137
+                end 
138
+            end
139
+            
140
+                reflec=reflec+RL;
141
+                reflec2=reflec2+(RL)^2;
142
+                transm=transm+TR;
143
+                transm2=transm2+(TR)^2;
144
+                for i=1:n+1
145
+                    SF(i)=SF(i)+SCAL(i);
146
+                    SF2(i)=SF2(i)+(SCAL(i))^2;
147
+                end
148
+            
149
+            
150
+            a=a+1;
151
+          
152
+end
153
+a=a-1;
154
+SF=SF./a;
155
+SF2=SF2./a;
156
+dlmwrite('SF_m.txt',SF,'delimiter',',','precision',18);
157
+kk=T/n;
158
+p=0:kk:T;                   % p is line along x-axis.
159
+plot(p,SF,'r'); hold on
160
+plot(p,SF2,'b');
161
+
162
+save per_bench.mat
163
+

+ 382
- 0
Vanilla MCode/New_SN_bench_solver_MC.m View File

@@ -0,0 +1,382 @@
1
+%THIS IS THE AUXILIARY FUNCTION OF THE MAIN CODE "New_SN_bench"
2
+
3
+%RICHARD VASQUES & NITIN KUMAR YADAV
4
+
5
+
6
+function[Z,extra,n1,randseed] = New_SN_bench_solver_MC(T,m1,m2,n,N,Es1,Es2,Et1,Et2,yo,y_,Q1,Q2,u,wt,randseed,med,a)
7
+
8
+%Building realization----------
9
+%EXP. Random Medium....
10
+if med==0
11
+    rng(randseed);
12
+    m12=zeros(2,1);
13
+    m12(1)=m1;m12(2)=m2;
14
+    xx=rand(1);
15
+    if xx<=m1/(m1+m2)
16
+        mm=0;
17
+    else
18
+        mm=1;
19
+    end
20
+    s=0;i=0;
21
+    while s<T
22
+        i=i+1;
23
+        x1=exprnd(m12(mod(mm,2)+1),1,1);
24
+        s=s+x1;
25
+        if s<=T
26
+            x(i,1)=s;
27
+        else
28
+            x(i,1)=T;
29
+        end
30
+        x(i,2)=mod(mm,2)+1;
31
+        mm=mm+1;
32
+    end
33
+    randseed=rng;
34
+else
35
+    interval=T/n;
36
+    %Periodic Medium........
37
+    m12=zeros(2,1);
38
+    m12(1)=m1;m12(2)=m2;
39
+    mm=0;
40
+    i=0;
41
+    s=0;
42
+    if a>1
43
+	if a<m1/interval+2
44
+		i=i+1;
45
+		x1=(a-1)*interval;
46
+		s=s+x1;
47
+		x(i,1)=s;
48
+		x(i,2)=2;
49
+	else
50
+    		i=i+1;
51
+		x1=(a-m1/interval-1)*interval;
52
+		s=s+x1;
53
+		x(i,1)=s;
54
+		x(i,2)=mod(mm,2)+1;
55
+        	mm=mm+1;
56
+	end
57
+    end
58
+    while s<T
59
+        i=i+1;
60
+        x1=m12(mod(mm,2)+1);
61
+        s=s+x1;
62
+        if s<=T
63
+            x(i,1)=s;
64
+        else
65
+            x(i,1)=T;
66
+        end
67
+        x(i,2)=mod(mm,2)+1;
68
+        mm=mm+1;
69
+    end
70
+end
71
+%Adding interfaces and extra points inside layers--------------
72
+
73
+%L(*,1) is the spatial point
74
+%L(*,2) indicates in which material is the spatial point (1 or 2).
75
+%If L(*,2)=3, the point is an interface.
76
+H=T/n;n1=1;i=1;j=1;
77
+extra=zeros(n,1);
78
+t1=i*H;
79
+L(1,1)=0;L(1,2)=x(j,2);
80
+if t1>x(j,1)
81
+    extra(i)=2;
82
+    h(n1)=x(j,1)/2;
83
+    L(n1+1,1)=h(n1);
84
+    L(n1+1,2)=x(j,2);
85
+    h(n1+1)=h(n1);
86
+    L(n1+2,1)=x(j,1);
87
+    L(n1+2,2)=3;
88
+    n1=n1+2;
89
+else
90
+    while t1<=x(j,1)
91
+        h(n1)=H;
92
+        L(n1+1,1)=i*H;
93
+        if L(n1+1,1)==x(j,1)                 % checking is point is interface point?
94
+            L(n1+1,2)=3;
95
+        else
96
+            L(n1+1,2)=x(j,2);
97
+        end
98
+        n1=n1+1;
99
+        i=i+1;
100
+        t1=i*H;
101
+    end
102
+    if x(j,1)<T && L(n1)~=x(j,1)
103
+        extra(i)=1;
104
+        L(n1+1,1)=x(j,1);
105
+        L(n1+1,2)=3;
106
+        h(n1)=L(n1+1,1)-L(n1,1);
107
+        n1=n1+1;
108
+    end
109
+    if t1==T && (x(j+1,1)-x(j,1))<H
110
+        extra(i)=extra(i)+1;
111
+        j=j+1;
112
+        L(n1+1,1)=(x(j,1)+x(j-1,1))/2;
113
+        L(n1+1,2)=x(j-1,2);
114
+        h(n1)=L(n1+1,1)-L(n1,1);
115
+        n1=n1+1;
116
+        L(n1+1,1)=x(j,1);
117
+        L(n1+1,2)=3;
118
+        h(n1)=L(n1+1,1)-L(n1,1);
119
+        n1=n1+1;
120
+        i=i+1;
121
+    end
122
+end
123
+j=2;
124
+if x(1,1)~=T
125
+    while i<=n
126
+        if t1>x(j,1)
127
+            extra(i)=extra(i)+2;
128
+            L(n1+1,1)=(x(j,1)+x(j-1,1))/2;
129
+            L(n1+1,2)=x(j,2);
130
+            h(n1)=L(n1+1,1)-L(n1,1);
131
+            n1=n1+1;
132
+            L(n1+1,1)=x(j,1);
133
+            L(n1+1,2)=3;
134
+            h(n1)=L(n1+1,1)-L(n1,1);
135
+            n1=n1+1;
136
+        else
137
+            while t1<=x(j,1)
138
+                L(n1+1,1)=i*H;
139
+                if L(n1+1,1)==x(j,1)          %checking is point is interface point?
140
+                    L(n1+1,2)=3;
141
+                else
142
+                    L(n1+1,2)=x(j,2);
143
+                end
144
+                h(n1)=L(n1+1,1)-L(n1,1);
145
+                n1=n1+1;
146
+                i=i+1;
147
+                t1=i*H;
148
+            end
149
+            if x(j,1)<T && L(n1)~=x(j,1)
150
+                extra(i)=extra(i)+1;
151
+                L(n1+1,1)=x(j,1);
152
+                L(n1+1,2)=3;
153
+                h(n1)=L(n1+1,1)-L(n1,1);
154
+                n1=n1+1;
155
+            end
156
+        end
157
+        j=j+1;
158
+        
159
+        if t1==T && (x(j,1)-x(j-1,1))<H
160
+            extra(i)=extra(i)+1;
161
+            L(n1+1,1)=(x(j,1)+x(j-1,1))/2;
162
+            L(n1+1,2)=x(j-1,2);
163
+            h(n1)=L(n1+1,1)-L(n1,1);
164
+            n1=n1+1;
165
+            L(n1+1,1)=x(j,1);
166
+            L(n1+1,2)=3;
167
+            h(n1)=L(n1+1,1)-L(n1,1);
168
+            n1=n1+1;
169
+            i=i+1;
170
+        end
171
+    end
172
+end
173
+n1=n1-1;
174
+
175
+A=zeros(N*n1,N*n1); B=zeros(N*n1,1);
176
+Et12=zeros(2,1); Es12=zeros(2,1); Q12=zeros(2,1);
177
+Et12(1)=Et1;Et12(2)=Et2;
178
+Es12(1)=Es1;Es12(2)=Es2;
179
+Q12(1)=Q1;Q12(2)=Q2;
180
+minleft=T/2-m1/2;
181
+maxright=T/2+m1/2;
182
+
183
+% Diagonal Block of matrix up to N/2.................................
184
+for t=1:N/2
185
+    s=(t-1)*n1;
186
+    A(s+1,s+1)=-u(t)*(1/h(1)+1/(h(1)+h(2)))+Et12(L(1,2))-Es12(L(1,2))*wt(t)/2;
187
+    A(s+1,s+2)=u(t)*(1/h(1)+1/h(2));
188
+    A(s+1,s+3)=-u(t)*(h(1)/(h(2)*(h(1)+h(2))));
189
+    %B(s+1)=Q12(L(1,2))/2;
190
+    B(s+1)=0.;
191
+    for i=2:n1-1
192
+        if L(i,2)==3
193
+            if i==n1-1
194
+                A(s+i,s+i)=-u(t)*(1/h(i)+1/(h(i)+h(i+1)))+Et12(L(i+1,2))-Es12(L(i+1,2))*wt(t)/2;
195
+                A(s+i,s+i+1)=u(t)*(1/h(i)+1/h(i+1));
196
+		if L(i+1,1)>=minleft && L(i+1,1)<=maxright 
197
+                	B(s+i)=u(t)*y_*(h(i)/(h(i+1)*(h(i)+h(i+1))))+Q12(L(i+1,2))/2;
198
+		else
199
+			B(s+i)=u(t)*y_*(h(i)/(h(i+1)*(h(i)+h(i+1))));
200
+		end
201
+            else
202
+                A(s+i,s+i)=-u(t)*(1/h(i)+1/(h(i)+h(i+1)))+Et12(L(i+1,2))-Es12(L(i+1,2))*wt(t)/2;
203
+                A(s+i,s+i+1)=u(t)*(1/h(i)+1/h(i+1));
204
+                A(s+i,s+i+2)=-u(t)*(h(i)/(h(i+1)*(h(i)+h(i+1))));
205
+		if L(i+1,1)>=minleft && L(i+1,1)<=maxright
206
+                	B(s+i)=Q12(L(i+1,2))/2;
207
+		else
208
+			B(s+i)=0.;
209
+		end
210
+            end
211
+        else
212
+            A(s+i,s+i-1)=-u(t)*h(i)/(h(i-1)*(h(i-1)+h(i)));
213
+            A(s+i,s+i)=u(t)*(h(i)-h(i-1))/(h(i)*h(i-1))+Et12(L(i,2))-Es12(L(i,2))*wt(t)/2;
214
+            A(s+i,s+i+1)=u(t)*h(i-1)/(h(i)*(h(i-1)+h(i)));
215
+	    if L(i,1)>=minleft && L(i,1)<=maxright
216
+            	B(s+i)=Q12(L(i,2))/2;
217
+	    else
218
+		B(s+i)=0.;
219
+	    end
220
+        end
221
+    end
222
+    A(s+n1,s+n1-1)=-u(t)*h(n1)/(h(n1-1)*(h(n1-1)+h(n1)));
223
+    A(s+n1,s+n1)=u(t)*(h(n1)-h(n1-1))/(h(n1)*h(n1-1))+Et12(L(n1,2))-Es12(L(n1,2))*wt(t)/2;
224
+    %B(s+n1)=-u(t)*y_*h(n1-1)/(h(n1)*(h(n1-1)+h(n1)))+Q12(L(n1,2))/2;
225
+    B(s+n1)=-u(t)*y_*h(n1-1)/(h(n1)*(h(n1-1)+h(n1)));
226
+    
227
+    % Remaining Blocks in same direction up to N/2..............
228
+    l=t;
229
+    if l==1 && N>2
230
+        for p=l+1:N/2
231
+            S=(p-1)*n1;
232
+            for i=1:n1
233
+                if L(i,2)==3
234
+                    A(s+i,S+i)=-Es12(L(i+1,2))*wt(p)/2;
235
+                else
236
+                    A(s+i,S+i)=-Es12(L(i,2))*wt(p)/2;
237
+                end
238
+            end
239
+        end
240
+    elseif l>1 && N>2
241
+        for p=1:l-1
242
+            S=(p-1)*n1;
243
+            for i=1:n1
244
+                if L(i,2)==3
245
+                    A(s+i,S+i)=-Es12(L(i+1,2))*wt(p)/2;
246
+                else
247
+                    A(s+i,S+i)=-Es12(L(i,2))*wt(p)/2;
248
+                end
249
+            end
250
+        end
251
+        for p=l+1:N/2
252
+            S=(p-1)*n1;
253
+            for i=1:n1
254
+                if L(i,2)==3
255
+                    A(s+i,S+i)=-Es12(L(i+1,2))*wt(p)/2;
256
+                else
257
+                    A(s+i,S+i)=-Es12(L(i,2))*wt(p)/2;
258
+                end;
259
+            end
260
+        end
261
+      
262
+    end
263
+    
264
+    % Blocks from N/2 to N........................................
265
+    a=0;
266
+    for p=1:N/2
267
+        S=(N/2+p-1)*n1;
268
+        for i=2:n1
269
+            if L(i,2)==3
270
+                A(s+i,S+i-1)=-Es12(L(i+1,2))*wt(N/2+p)/2;
271
+            else
272
+                A(s+i,S+i-1)=-Es12(L(i,2))*wt(N/2+p)/2;
273
+            end
274
+        end
275
+        a=a+(Es12(L(1,2))*wt(N/2+p)*yo/2);
276
+    end
277
+    B(s+1)=B(s+1)+a;
278
+end
279
+    
280
+% Diagonal Block of matrix from N/2+1 to N.........................
281
+for t=N/2+1:N
282
+    s=(t-1)*n1;
283
+    A(s+1,s+1)=u(t)*(h(2)-h(1))/(h(2)*h(1))+Et12(L(1,2))-Es12(L(1,2))*wt(t)/2;
284
+    A(s+1,s+1+1) = u(t)*h(1)/(h(2)*(h(1)+h(2)));
285
+    %B(s+1)=u(t)*yo*h(2)/(h(1)*(h(1)+h(2)))+Q12(L(1,2))/2;
286
+    B(s+1)=u(t)*yo*h(2)/(h(1)*(h(1)+h(2)));
287
+    for i=2:n1-1
288
+        if L(i+1,2)==3
289
+            if i==2
290
+                A(s+i,s+i)=u(t)*(1/h(i)+1/(h(i)+h(i-1)))+Et12(L(i,2))-Es12(L(i,2))*wt(t)/2;
291
+                A(s+i,s+i-1)=-u(t)*(1/h(i)+1/h(i-1));
292
+		if L(i,1)>=minleft && L(i,1)<=maxright
293
+                	B(s+i)=-u(t)*yo*(h(i)/(h(i-1)*(h(i)+h(i-1))))+Q12(L(i,2))/2;
294
+		else
295
+			B(s+i)=-u(t)*yo*(h(i)/(h(i-1)*(h(i)+h(i-1))));
296
+		end
297
+            else
298
+                A(s+i,s+i)=u(t)*(1/h(i)+1/(h(i)+h(i-1)))+Et12(L(i,2))-Es12(L(i,2))*wt(t)/2;
299
+                A(s+i,s+i-1)=-u(t)*(1/h(i)+1/h(i-1));
300
+                A(s+i,s+i-2)=u(t)*(h(i)/(h(i-1)*(h(i)+h(i-1))));
301
+		if L(i,1)>=minleft && L(i,1)<=maxright
302
+                	B(s+i)=Q12(L(i,2))/2;
303
+		else
304
+			B(s+i)=0.;
305
+		end
306
+            end
307
+        else
308
+            A(s+i,s+i-1)=-u(t)*h(i+1)/(h(i)*(h(i)+h(i+1)));
309
+            A(s+i,s+i)=u(t)*(h(i+1)-h(i))/(h(i+1)*h(i))+Et12(L(i+1,2))-Es12(L(i+1,2))*wt(t)/2;
310
+            A(s+i,s+i+1)=u(t)*h(i)/(h(i+1)*(h(i)+h(i+1)));
311
+	    if L(i+1,1)>=minleft && L(i+1,1)<=maxright
312
+            	B(s+i)=Q12(L(i+1,2))/2;
313
+	    else
314
+		B(s+i)=0.;
315
+	    end
316
+        end
317
+    end
318
+    A(s+n1,s+n1)=u(t)*(1/h(n1)+1/(h(n1)+h(n1-1)))+Et12(L(n1,2))-Es12(L(n1,2))*wt(t)/2;
319
+    A(s+n1,s+n1-1)=-u(t)*(1/h(n1)+1/h(n1-1));
320
+    A(s+n1,s+n1-2)=u(t)*(h(n1)/(h(n1-1)*(h(n1)+h(n1-1))));
321
+    %B(s+n1)=Q12(L(n1,2))/2;
322
+    B(s+n1)=0.;
323
+    
324
+    % Remaining Blocks in same direction up to N..............
325
+    l=t;
326
+    if l==N/2+1 && N>2
327
+        for p=l+1:N
328
+            S=(p-1)*n1;
329
+            for i=1:n1
330
+                if L(i+1,2)==3
331
+                    A(s+i,S+i)=-Es12(L(i,2))*wt(p)/2;
332
+                else
333
+                    A(s+i,S+i)=-Es12(L(i+1,2))*wt(p)/2;
334
+                end
335
+            end
336
+        end
337
+    elseif l>N/2+1 && N>2
338
+        for p=N/2+1:l-1
339
+            S=(p-1)*n1;
340
+            for i=1:n1
341
+                if L(i+1,2)==3
342
+                    A(s+i,S+i)=-Es12(L(i,2))*wt(p)/2;
343
+                else
344
+                    A(s+i,S+i)=-Es12(L(i+1,2))*wt(p)/2;
345
+                end
346
+            end
347
+        end
348
+        for p=l+1:N
349
+            S=(p-1)*n1;
350
+            for i=1:n1
351
+                if L(i+1,2)==3
352
+                    A(s+i,S+i)=-Es12(L(i,2))*wt(p)/2;
353
+                else
354
+                    A(s+i,S+i)=-Es12(L(i+1,2))*wt(p)/2;
355
+                end
356
+            end
357
+        end
358
+      
359
+    end
360
+    
361
+    % Blocks from 1 to N/2........................................
362
+    a=0;
363
+    for p=1:N/2
364
+        S=(p-1)*n1;
365
+        for i=1:n1-1
366
+            if L(i+1,2)==3
367
+                A(s+i,S+i+1)=-Es12(L(i,2))*wt(p)/2;
368
+            else
369
+                A(s+i,S+i+1)=-Es12(L(i+1,2))*wt(p)/2;
370
+            end
371
+        end
372
+        a=a+(Es12(L(n1,2))*wt(p)*y_/2);
373
+    end
374
+    B(s+n1)=B(s+n1)+a;
375
+end
376
+
377
+Z=A\B; % Solving for angular flux.
378
+return
379
+end
380
+
381
+%If you want to see Flux outgoing in positive directin see Y(n) or Z(M+1)
382
+%If you want see Flux scattering back see Y(n+1) or Z(M+2)

+ 212
- 0
Vanilla MCode/SN_hom.m View File

@@ -0,0 +1,212 @@
1
+%THIS CODE SOLVES THE STEADY-STATE, MONOENERGETIC TRANSPORT EQUATION 
2
+%IN A HOMOGENEOUS MEDIUM WITH ISOTROPIC SCATTERING AND ISOTROPIC SOURCE
3
+%IN ROD AND SLAB (SN) GEOMETRIES. IT USES 2-POINT CENTRAL DIFFERENCING
4
+%(ORDER H^2) WITH 3-POINT FORWARD/BACKWARD AT THE BOUNDARIES.
5
+
6
+%MAIN OUTPUTS: transmission "TR"; reflection "RL"
7
+%              scalar flux "SCAL"
8
+
9
+%RICHARD VASQUES & NITIN KUMAR YADAV
10
+
11
+clear all
12
+clc
13
+% INPUTS--------------------------------
14
+display('Default is ROD geometry!!');
15
+display('  ');
16
+rod_slab=input('Enter 1 if you want to change to SLAB Geometry:   ');
17
+if rod_slab==1
18
+    N=1;
19
+    while floor(N/2)~=ceil(N/2) 
20
+        N=input(   'Enter the number of Discrete ordinate directions: ');
21
+    end
22
+    clc
23
+    display('SOLVING PROBLEM IN SLAB GEOMETRY!!');
24
+    display('  ')
25
+else
26
+    rod_slab=0;
27
+    clc
28
+    display('SOLVING PROBLEM IN ROD GEOMETRY!!');
29
+    display('  ');
30
+    N=2;
31
+end
32
+T=input('Enter the total length of the system:               ');
33
+n=1;
34
+while n<4
35
+    n=input('Enter at how many points you want to calculate:     ');
36
+end
37
+Et=input('Enter the total cross section Sigma_t :             ');
38
+cc=10;
39
+while cc>1 || cc<0
40
+    cc=input('Enter the Scattering Ratio c (between 0 and 1):     ');
41
+end
42
+Es = cc*Et;
43
+yo=input('Enter the boundary value in the positive direction: ');
44
+y_=input('Enter the boundary value in the negative direction: ');
45
+Q=input('Enter the homogeneous isotropic Source:             ');
46
+
47
+M=n*N;  h=T/n;
48
+A=zeros(M,M);      B(1:M)=Q/2;
49
+
50
+% GAUSS-LEGENDRE QUADRATURE............................
51
+beta = (1:N-1)./sqrt(4*((1:N-1)).^2-1);
52
+[w,x] = eig(diag(beta,-1)+diag(beta,1));
53
+u = diag(x);
54
+wt = 2*w(1,:)'.^2;
55
+if rod_slab~=1
56
+    u(1)=-1;u(2)=1;
57
+end
58
+
59
+% Diagonal Block of matrix up to N/2.................................
60
+for t=1:N/2
61
+    s=(t-1)*n;
62
+    A(s+1,s+1)=(-11*u(t)/(6*h)+Et-Es*wt(t)/2);
63
+    A(s+1,s+2) = 3*u(t)/h;
64
+    A(s+1,s+3) = -3*u(t)/(2*h);
65
+    A(s+1,s+4) = u(t)/(3*h);
66
+    for i=2:n-1
67
+        A(s+i,s+i-1) = -u(t)/(2*h); 
68
+        A(s+i,s+i)=(Et-Es*wt(t)/2);
69
+        A(s+i,s+i+1)= u(t)/(2*h);
70
+    end
71
+    A(s+n,s+n-1) = -u(t)/(2*h);
72
+    A(s+n,s+n)=(Et-Es*wt(t)/2);
73
+    B(s+n)=-u(t)*y_/(2*h)+Q/2;
74
+    % Remaining Blocks in same direction up to N/2..............
75
+    l=t;
76
+    if l==1 && N>2
77
+        for p=l+1:N/2
78
+            S=(p-1)*n;
79
+            for i=1:n
80
+                A(s+i,S+i)=-Es*wt(p)/2;
81
+            end
82
+        end
83
+    elseif l>1 && N>2
84
+        for p=1:l-1
85
+            S=(p-1)*n;
86
+            for i=1:n
87
+                A(s+i,S+i)=-Es*wt(p)/2;
88
+            end
89
+        end
90
+        for p=l+1:N/2
91
+            S=(p-1)*n;
92
+            for i=1:n
93
+                A(s+i,S+i)=-Es*wt(p)/2;
94
+            end
95
+        end
96
+      
97
+    end
98
+    % Blocks from N/2 to N........................................
99
+    a=0;
100
+    for p=1:N/2
101
+        S=(N/2+p-1)*n;
102
+        for i=2:n
103
+            A(s+i,S+i-1)=-Es*wt(N/2+p)/2;
104
+        end
105
+        a=a+(Es*wt(N/2+p)*yo/2);
106
+    end
107
+    B(s+1)=a+Q/2;
108
+end
109
+
110
+
111
+% Diagonal Block of matrix from N/2+1 to N.........................
112
+for t=N/2+1:N
113
+    s=(t-1)*n;
114
+    A(s+1,s+1)=(Et-Es*wt(t)/2);
115
+    A(s+1,s+1+1) = u(t)/(2*h);
116
+    B(s+1)=u(t)*yo/(2*h)+Q/2;
117
+    for i=2:n-1
118
+        A(s+i,s+i-1) = -u(t)/(2*h); 
119
+        A(s+i,s+i)=(Et-Es*wt(t)/2);
120
+        A(s+i,s+i+1)= u(t)/(2*h);
121
+    end
122
+    A(s+n,s+n-3)=-u(t)/(3*h);
123
+    A(s+n,s+n-2)=3*u(t)/(2*h);
124
+    A(s+n,s+n-1)=-3*u(t)/h;
125
+    A(s+n,s+n) =(11*u(t)/(6*h)+Et-Es*wt(t)/2);
126
+    
127
+    % Remaining Blocks in same direction up to N..............
128
+    l=t;
129
+    if l==N/2+1 && N>2
130
+        for p=l+1:N
131
+            S=(p-1)*n;
132
+            for i=1:n
133
+                A(s+i,S+i)=-Es*wt(p)/2;
134
+            end
135
+        end
136
+    elseif l>N/2+1 && N>2
137
+        for p=N/2+1:l-1
138
+            S=(p-1)*n;
139
+            for i=1:n
140
+                A(s+i,S+i)=-Es*wt(p)/2;
141
+            end
142
+        end
143
+        for p=l+1:N
144
+            S=(p-1)*n;
145
+            for i=1:n
146
+                A(s+i,S+i)=-Es*wt(p)/2;
147
+            end
148
+        end
149
+      
150
+    end
151
+    
152
+    % Blocks from 1 to N/2........................................
153
+    a=0;
154
+    for p=1:N/2
155
+        S=(p-1)*n;
156
+        for i=1:n-1
157
+            A(s+i,S+i+1)=-Es*wt(p)/2;
158
+        end
159
+        a=a+(Es*wt(p)*y_/2);
160
+    end
161
+    B(s+n)=a+Q/2;
162
+end
163
+
164
+X=A\B';
165
+
166
+Y=zeros(M+N,1);
167
+
168
+% Adding boundary conditions to the array...............
169
+i=1;
170
+for t=1:N/2
171
+    s=(t-1);
172
+    for j=i:t*n+s
173
+        Y(j)=X(j-s);
174
+    end
175
+    Y(j+1)=y_;
176
+    i=j+2;
177
+end
178
+i=i-1;
179
+for t=N/2+1:N
180
+    Y(i+1)=yo;
181
+    i=i+1;
182
+    for j=i+1:t*n+t
183
+        Y(j)=X(j-t);
184
+    end 
185
+    i=j;
186
+end
187
+
188
+% Calculating Reflection, Transmission and Scalar Flux...........
189
+RL=0;
190
+for t=1:N/2
191
+    s=(t-1)*n;
192
+    RL=RL+abs(wt(t)*u(t)*Y(s+t));
193
+end
194
+TR=0;
195
+for t=N/2+1:N
196
+ s=(t-1)*n;
197
+ S=(t-N/2-1)*n;
198
+ k=N/2+1;
199
+    TR=TR+wt(t)*u(t)*Y(s+n+t);
200
+end
201
+m=n+1;
202
+SCAL=zeros(m,1);
203
+for t=1:N
204
+    s=(t-1)*m;
205
+    for i=1:m
206
+        SCAL(i)=SCAL(i)+wt(t)*Y(s+i);
207
+    end 
208
+end
209
+
210
+x=0:h:T;
211
+%plot(x,SCAL,'--');
212
+plot(x,SCAL,'b'); hold on

+ 207
- 0
Vanilla MCode/SN_per_bench_solver.m View File

@@ -0,0 +1,207 @@
1
+%THIS IS THE AUXILIARY FUNCTION OF THE MAIN CODE "periodic_benchmark"
2
+
3
+%RICHARD VASQUES & NITIN KUMAR YADAV
4
+
5
+
6
+function[Z,extra,n1,B] = SN_per_bench_solver(T,m1,m2,n,N,Es1,Es2,Et1,Et2,yo,y_,Q1,Q2,u,wt,a)
7
+
8
+%Periodic Medium........
9
+interval=T/n;
10
+m12=zeros(2,1);
11
+m12(1)=m1;m12(2)=m2;
12
+mm=0;
13
+i=0;
14
+s=0;
15
+if a>1
16
+    if a<m2/interval+2
17
+		i=i+1;
18
+		x1=(a-1)*interval;
19
+		s=s+x1;
20
+		x(i,1)=s;
21
+		x(i,2)=2;
22
+	else
23
+    	i=i+1;
24
+		x1=(a-m2/interval-1)*interval;
25
+		s=s+x1;
26
+		x(i,1)=s;
27
+		x(i,2)=mod(mm,2)+1;
28
+       mm=mm+1;
29
+    end
30
+end
31
+
32
+while s<T
33
+    i=i+1;
34
+    x1=m12(mod(mm,2)+1);
35
+    s=s+x1;
36
+    if s<=T
37
+        x(i,1)=s;
38
+    else
39
+        x(i,1)=T;
40
+    end
41
+    x(i,2)=mod(mm,2)+1;
42
+    mm=mm+1;
43
+end
44
+
45
+H=T/n;n1=1;i=1;j=1;
46
+extra=zeros(n,1);
47
+t1=i*H;
48
+L(1,1)=0;L(1,2)=x(j,2);
49
+if t1==x(j,1)
50
+    extra(i)=1;
51
+    h(n1)=x(j,1)/2;
52
+    L(n1+1,1)=h(n1);
53
+    L(n1+1,2)=x(j,2);
54
+    h(n1+1)=h(n1);
55
+    L(n1+2,1)=x(j,1);
56
+    L(n1+2,2)=3;
57
+    n1=n1+2;
58
+     i=i+1;
59
+     t1=i*H;
60
+else
61
+    while t1<=x(j,1)
62
+        h(n1)=H;
63
+        L(n1+1,1)=i*H;
64
+        if L(n1+1,1)==x(j,1)                 % checking is point is interface point?
65
+            L(n1+1,2)=3;
66
+        else
67
+            L(n1+1,2)=x(j,2);
68
+        end
69
+        n1=n1+1;
70
+        i=i+1;
71
+        t1=i*H;
72
+    end
73
+end
74
+j=2;
75
+if x(1,1)~=T
76
+    while i<=n
77
+        while t1<=x(j,1)
78
+                L(n1+1,1)=i*H;
79
+                if L(n1+1,1)==x(j,1)          %checking is point is interface point?
80
+                    L(n1+1,2)=3;
81
+                else
82
+                    L(n1+1,2)=x(j,2);
83
+                end
84
+                h(n1)=L(n1+1,1)-L(n1,1);
85
+                n1=n1+1;
86
+                i=i+1;
87
+                t1=i*H;
88
+        end
89
+        j=j+1;
90
+        if L(n1,1)==T && L(n1,2)==3 && L(n1-1,2)==3
91
+            i=i-1; 
92
+            extra(i)=extra(i)+1;
93
+            L(n1,1)=(x(j-1,1)+x(j-2,1))/2;
94
+            L(n1,2)=x(j-1,2);
95
+            h(n1-1)=L(n1,1)-L(n1-1,1);
96
+            n1=n1+1;
97
+            L(n1,1)=x(j-1,1);
98
+            L(n1,2)=3;
99
+            h(n1-1)=L(n1,1)-L(n1-1,1);
100
+            i=i+1;
101
+        end       
102
+    end
103
+end
104
+n1=n1-1;
105
+
106
+A=zeros(N*n1,N*n1); B=zeros(N*n1,1);
107
+Et12=zeros(2,1); Es12=zeros(2,1); Q12=zeros(2,1);
108
+Et12(1)=Et1;Et12(2)=Et2;
109
+Es12(1)=Es1;Es12(2)=Es2;
110
+Q12(1)=Q1;Q12(2)=Q2;
111
+
112
+
113
+% Diagonal Block of matrix up to N/2.................................
114
+for t=1:N/2
115
+    s=(t-1)*n1;
116
+    A(s+1,s+1)=-u(t)*(1/h(1))+Et12(L(1,2))-Es12(L(1,2))*wt(t)/2;
117
+    A(s+1,s+2)=u(t)*(1/h(1));
118
+    B(s+1)=Q12(L(1,2))/2;
119
+    for i=2:n1-1
120
+        if L(i,2)==3
121
+            if i==n1-1
122
+                A(s+i,s+i)=-u(t)*(1/h(i)+1/(h(i)+h(i+1)))+Et12(L(i+1,2))-Es12(L(i+1,2))*wt(t)/2;
123
+                A(s+i,s+i+1)=u(t)*(1/h(i)+1/h(i+1));
124
+                B(s+i)=u(t)*y_*(h(i)/(h(i+1)*(h(i)+h(i+1))))+Q12(L(i+1,2))/2;
125
+            else
126
+                A(s+i,s+i)=-u(t)*(1/h(i)+1/(h(i)+h(i+1)))+Et12(L(i+1,2))-Es12(L(i+1,2))*wt(t)/2;
127
+                A(s+i,s+i+1)=u(t)*(1/h(i)+1/h(i+1));
128
+                A(s+i,s+i+2)=-u(t)*(h(i)/(h(i+1)*(h(i)+h(i+1))));
129
+                B(s+i)=Q12(L(i+1,2))/2;
130
+        
131
+            end
132
+        else
133
+            A(s+i,s+i-1)=-u(t)*h(i)/(h(i-1)*(h(i-1)+h(i)));
134
+            A(s+i,s+i)=u(t)*(h(i)-h(i-1))/(h(i)*h(i-1))+Et12(L(i,2))-Es12(L(i,2))*wt(t)/2;
135
+            A(s+i,s+i+1)=u(t)*h(i-1)/(h(i)*(h(i-1)+h(i)));
136
+            	B(s+i)=Q12(L(i,2))/2;
137
+        end
138
+    end
139
+    A(s+n1,s+n1-1)=-u(t)*h(n1)/(h(n1-1)*(h(n1-1)+h(n1)));
140
+    A(s+n1,s+n1)=u(t)*(h(n1)-h(n1-1))/(h(n1)*h(n1-1))+Et12(L(n1,2))-Es12(L(n1,2))*wt(t)/2;
141
+    B(s+n1)=-u(t)*y_*h(n1-1)/(h(n1)*(h(n1-1)+h(n1)))+Q12(L(n1,2))/2;
142
+    
143
+    % Blocks from N/2 to N........................................
144
+    a=0;
145
+    for p=1:N/2
146
+        S=(N/2+p-1)*n1;
147
+        for i=2:n1
148
+            if L(i,2)==3
149
+                A(s+i,S+i-1)=-Es12(L(i+1,2))*wt(N/2+p)/2;
150
+            else
151
+                A(s+i,S+i-1)=-Es12(L(i,2))*wt(N/2+p)/2;
152
+            end
153
+        end
154
+        a=a+(Es12(L(1,2))*wt(N/2+p)*yo/2);
155
+    end
156
+    B(s+1)=B(s+1)+a;
157
+end
158
+    
159
+% Diagonal Block of matrix from N/2+1 to N.........................
160
+for t=N/2+1:N
161
+    s=(t-1)*n1;
162
+    A(s+1,s+1)=u(t)*(h(2)-h(1))/(h(2)*h(1))+Et12(L(1,2))-Es12(L(1,2))*wt(t)/2;
163
+    A(s+1,s+1+1) = u(t)*h(1)/(h(2)*(h(1)+h(2)));
164
+    B(s+1)=u(t)*yo*h(2)/(h(1)*(h(1)+h(2)))+Q12(L(1,2))/2;
165
+    for i=2:n1-1
166
+        if L(i+1,2)==3
167
+            if i==2
168
+                A(s+i,s+i)=u(t)*(1/h(i)+1/(h(i)+h(i-1)))+Et12(L(i,2))-Es12(L(i,2))*wt(t)/2;
169
+                A(s+i,s+i-1)=-u(t)*(1/h(i)+1/h(i-1));
170
+                B(s+i)=-u(t)*yo*(h(i)/(h(i-1)*(h(i)+h(i-1))))+Q12(L(i,2))/2;
171
+            else
172
+                A(s+i,s+i)=u(t)*(1/h(i)+1/(h(i)+h(i-1)))+Et12(L(i,2))-Es12(L(i,2))*wt(t)/2;
173
+                A(s+i,s+i-1)=-u(t)*(1/h(i)+1/h(i-1));
174
+                A(s+i,s+i-2)=u(t)*(h(i)/(h(i-1)*(h(i)+h(i-1))));
175
+                B(s+i)=Q12(L(i,2))/2;
176
+	        end
177
+        else
178
+            A(s+i,s+i-1)=-u(t)*h(i+1)/(h(i)*(h(i)+h(i+1)));
179
+            A(s+i,s+i)=u(t)*(h(i+1)-h(i))/(h(i+1)*h(i))+Et12(L(i+1,2))-Es12(L(i+1,2))*wt(t)/2;
180
+            A(s+i,s+i+1)=u(t)*h(i)/(h(i+1)*(h(i)+h(i+1)));
181
+            B(s+i)=Q12(L(i+1,2))/2;
182
+	    end
183
+    end
184
+    A(s+n1,s+n1)=u(t)*(1/h(n1))+Et12(L(n1,2))-Es12(L(n1,2))*wt(t)/2;
185
+    A(s+n1,s+n1-1)=-u(t)*(1/h(n1));
186
+    B(s+n1)=Q12(L(n1,2))/2;
187
+    
188
+    % Blocks from 1 to N/2........................................
189
+    a=0;
190
+    for p=1:N/2
191
+        S=(p-1)*n1;
192
+        for i=1:n1-1
193
+            if L(i+1,2)==3
194
+                A(s+i,S+i+1)=-Es12(L(i,2))*wt(p)/2;
195
+            else
196
+                A(s+i,S+i+1)=-Es12(L(i+1,2))*wt(p)/2;
197
+            end
198
+        end
199
+        a=a+(Es12(L(n1,2))*wt(p)*y_/2);
200
+    end
201
+    B(s+n1)=B(s+n1)+a;
202
+end
203
+
204
+Z=A\B; % Solving for angular flux.
205
+
206
+return
207
+end

+ 143
- 0
Vanilla MCode/periodic_benchmark.m View File

@@ -0,0 +1,143 @@
1
+%THIS CODE GENERATES ENSEMBLE-AVERAGED BENCHMARK RESULTS FOR THE
2
+%STEADY-STATE, MONOENERGETIC TRANSPORT EQUATION IN A PERIODIC MEDIUM
3
+%COMPOSED OF TWO MATERIALS WITH ISOTROPIC SCATTERING AND ISOTROPIC SOURCE
4
+%IN ROD GEOMETRY. THE TOTAL NUMBER OF REALIZATIONS IS DEFINED BY THE
5
+%SPATIAL DISCRETIZATION.
6
+
7
+%THIS CODE CALLS THE SOLVER ROUTINE "SN_per_bench_solver", USING 2-POINT CENTRAL
8
+%DIFFERENCING (ORDER H^2).
9
+
10
+clear all
11
+clc
12
+
13
+%------------------------------------------------------------
14
+%INPUTS
15
+T= 20;       % total length of the system
16
+n= 2560;     % # of points (discretization)
17
+
18
+yo=0.;      % left boundary condition
19
+y_=0.;      % right boundary condition
20
+
21
+m1= 1;    % thickness of material 1 layers
22
+m2= 1;    % thickness of material 2 layers
23
+
24
+Et1= 1.0;    % total cross section Sigma_t1 of material 1
25
+cc= 0.5;    % Scattering Ratio c1 of material 1
26
+Q1= 1.0;     % homogeneous isotropic Source of material 1
27
+
28
+Et2=0.;      % total cross section Sigma_t2 of material 2
29
+cc2= 0.;    % Scattering Ratio c2 of material 2
30
+Q2= 0.;     % homogeneous isotropic Source of material 2
31
+%------------------------------------------------------------
32
+%------------------------------------------------------------
33
+
34
+Es1=cc*Et1;
35
+Ea1= Et1-Es1;
36
+Es2=cc2*Et2;
37
+Ea2= Et2-Es2;
38
+
39
+%weights and directions
40
+N=2;
41
+wt(1)= 1;
42
+wt(2)=1;
43
+u(1)=-1;
44
+u(2)=1;
45
+
46
+a=1;
47
+reflec=0;reflec2=0;
48
+transm=0;transm2=0;
49
+SF=zeros(n+1,1);
50
+SF2=zeros(n+1,1);
51
+cond=0;
52
+
53
+%MAIN LOOP
54
+total = (m1+m2)/(T/n);
55
+while(a<total+1)
56
+    fprintf('Problem %d of %d\n', a, total);
57
+    [Z,extra,n1,B] = SN_per_bench_solver(T,m1,m2,n,N,Es1,Es2,Et1,Et2,yo,y_,Q1,Q2,u,wt,a);
58
+    
59
+            % Adjusting points..........................
60
+            X=zeros(n*N,1);
61
+            for i=1:N/2
62
+                X((i-1)*n+1)=Z((i-1)*n1+1);
63
+                k=2;
64
+                for j=2:n
65
+                    k=k+extra(j-1);
66
+                    X((i-1)*n+j)=Z((i-1)*n1+k);
67
+                    k=k+1;
68
+                end
69
+            end            
70
+            for i=N/2+1:N
71
+                k=1;
72
+                for j=1:n
73
+                    k=k+extra(j);
74
+                    X((i-1)*n+j)=Z((i-1)*n1+k);
75
+                    k=k+1;
76
+                end
77
+            end
78
+            
79
+            Y=zeros(N*(n+1),1);
80
+            % Adding boundary conditions...............
81
+            i=1;
82
+            for t=1:N/2
83
+                s=(t-1);
84
+                for j=i:t*n+s
85
+                    Y(j)=X(j-s);
86
+                end
87
+                Y(j+1)=y_;
88
+                i=j+2;
89
+            end
90
+            i=i-1;
91
+            for t=N/2+1:N
92
+                Y(i+1)=yo;
93
+                i=i+1;
94
+                for j=i+1:t*n+t
95
+                    Y(j)=X(j-t);
96
+                end 
97
+                i=j;
98
+            end
99
+ 
100
+            % Calculating Reflection, Transmission and Scalar Flux......     
101
+            RL=0;
102
+            for t=1:N/2
103
+                s=(t-1)*n;
104
+                RL=RL+abs(wt(t)*u(t)*Y(s+t));
105
+            end
106
+            TR=0;
107
+            for t=N/2+1:N
108
+                s=(t-1)*n;
109
+                TR=TR+wt(t)*u(t)*Y(s+n+t);
110
+            end
111
+            m=n+1;
112
+            SCAL=zeros(m,1);
113
+            for t=1:N
114
+                s=(t-1)*m;
115
+                for i=1:m
116
+                    SCAL(i)=SCAL(i)+wt(t)*Y(s+i);
117
+                end 
118
+            end
119
+            
120
+                reflec=reflec+RL;
121
+                reflec2=reflec2+(RL)^2;
122
+                transm=transm+TR;
123
+                transm2=transm2+(TR)^2;
124
+                for i=1:n+1
125
+                    SF(i)=SF(i)+SCAL(i);
126
+                    SF2(i)=SF2(i)+(SCAL(i))^2;
127
+                end
128
+            
129
+            
130
+            a=a+1;
131
+          
132
+end
133
+a=a-1;
134
+SF=SF./a;
135
+SF2=SF2./a;
136
+
137
+kk=T/n;
138
+p=0:kk:T;                   % p is line along x-axis.
139
+plot(p,SF,'r'); hold on
140
+plot(p,SF2,'b');
141
+
142
+save per_bench.mat
143
+

+ 32
- 0
compare.py View File

@@ -0,0 +1,32 @@
1
+import csv
2
+import itertools
3
+
4
+var = "extra"
5
+
6
+file1 = open(var + ".csv", "rb")
7
+file2 = open("/Users/clayshieh/Documents/MATLAB/" + var + "_m.txt", "rb")
8
+
9
+r1 = csv.reader(file1)
10
+r2 = csv.reader(file2)
11
+row_count1 = sum(1 for row in r1)
12
+row_count2 = sum(1 for row in r2)
13
+if row_count1 != row_count2:
14
+	print "Actual has " + str(row_count1) + " rows Expected has " + str(row_count2) + "rows."
15
+file1.close()
16
+file2.close()
17
+file1 = open(var + ".csv", "rb")
18
+file2 = open("/Users/clayshieh/Documents/MATLAB/" + var + "_m.txt", "rb")
19
+reader1 = csv.reader(file1)
20
+reader2 = csv.reader(file2)
21
+
22
+
23
+row = 0
24
+for lhs, rhs in zip(reader1, reader2):
25
+	column = 0
26
+	for lh, rh in zip(lhs, rhs):
27
+		if abs(float(lh) - float(rh)) >= 1e-10:
28
+			print "Actual: " + str(lh) + " Expected: " + str(rh) + "   at: (" + str(row) + "," + str(column) + ")"
29
+		column += 1
30
+
31
+	row += 1
32
+print row

+ 202
- 0
csv/A.csv
File diff suppressed because it is too large
View File


+ 202
- 0
csv/B.csv View File

@@ -0,0 +1,202 @@
1
+5.000000000000000000e-01
2
+5.000000000000000000e-01
3
+5.000000000000000000e-01
4
+5.000000000000000000e-01
5
+0.000000000000000000e+00
6
+0.000000000000000000e+00
7
+0.000000000000000000e+00
8
+0.000000000000000000e+00
9
+0.000000000000000000e+00
10
+5.000000000000000000e-01
11
+5.000000000000000000e-01
12
+5.000000000000000000e-01
13
+5.000000000000000000e-01
14
+5.000000000000000000e-01
15
+0.000000000000000000e+00
16
+0.000000000000000000e+00
17
+0.000000000000000000e+00
18
+0.000000000000000000e+00
19
+0.000000000000000000e+00
20
+5.000000000000000000e-01
21
+5.000000000000000000e-01
22
+5.000000000000000000e-01
23
+5.000000000000000000e-01
24
+5.000000000000000000e-01
25
+0.000000000000000000e+00
26
+0.000000000000000000e+00
27
+0.000000000000000000e+00
28
+0.000000000000000000e+00
29
+0.000000000000000000e+00
30
+5.000000000000000000e-01
31
+5.000000000000000000e-01
32
+5.000000000000000000e-01
33
+5.000000000000000000e-01
34
+5.000000000000000000e-01
35
+0.000000000000000000e+00
36
+0.000000000000000000e+00
37
+0.000000000000000000e+00
38
+0.000000000000000000e+00
39
+0.000000000000000000e+00
40
+5.000000000000000000e-01
41
+5.000000000000000000e-01
42
+5.000000000000000000e-01
43
+5.000000000000000000e-01
44
+5.000000000000000000e-01
45
+0.000000000000000000e+00
46
+0.000000000000000000e+00
47
+0.000000000000000000e+00
48
+0.000000000000000000e+00
49
+0.000000000000000000e+00
50
+5.000000000000000000e-01
51
+5.000000000000000000e-01
52
+5.000000000000000000e-01
53
+5.000000000000000000e-01
54
+5.000000000000000000e-01
55
+0.000000000000000000e+00
56
+0.000000000000000000e+00
57
+0.000000000000000000e+00
58
+0.000000000000000000e+00
59
+0.000000000000000000e+00
60
+5.000000000000000000e-01
61
+5.000000000000000000e-01
62
+5.000000000000000000e-01
63
+5.000000000000000000e-01
64
+5.000000000000000000e-01
65
+0.000000000000000000e+00
66
+0.000000000000000000e+00
67
+0.000000000000000000e+00
68
+0.000000000000000000e+00
69
+0.000000000000000000e+00
70
+5.000000000000000000e-01
71
+5.000000000000000000e-01
72
+5.000000000000000000e-01
73
+5.000000000000000000e-01
74
+5.000000000000000000e-01
75
+0.000000000000000000e+00
76
+0.000000000000000000e+00
77
+0.000000000000000000e+00
78
+0.000000000000000000e+00
79
+0.000000000000000000e+00
80
+5.000000000000000000e-01
81
+5.000000000000000000e-01
82
+5.000000000000000000e-01
83
+5.000000000000000000e-01
84
+5.000000000000000000e-01
85
+0.000000000000000000e+00
86
+0.000000000000000000e+00
87
+0.000000000000000000e+00
88
+0.000000000000000000e+00
89
+0.000000000000000000e+00
90
+5.000000000000000000e-01
91
+5.000000000000000000e-01
92
+5.000000000000000000e-01
93
+5.000000000000000000e-01
94
+5.000000000000000000e-01
95
+0.000000000000000000e+00
96
+0.000000000000000000e+00
97
+0.000000000000000000e+00
98
+0.000000000000000000e+00
99
+0.000000000000000000e+00
100
+5.000000000000000000e-01
101
+5.000000000000000000e-01
102
+5.000000000000000000e-01
103
+5.000000000000000000e-01
104
+5.000000000000000000e-01
105
+5.000000000000000000e-01
106
+0.000000000000000000e+00
107
+0.000000000000000000e+00
108
+0.000000000000000000e+00
109
+0.000000000000000000e+00
110
+0.000000000000000000e+00
111
+5.000000000000000000e-01
112
+5.000000000000000000e-01
113
+5.000000000000000000e-01
114
+5.000000000000000000e-01
115
+0.000000000000000000e+00
116
+0.000000000000000000e+00
117
+0.000000000000000000e+00
118
+0.000000000000000000e+00
119
+0.000000000000000000e+00
120
+5.000000000000000000e-01
121
+5.000000000000000000e-01
122
+5.000000000000000000e-01
123
+5.000000000000000000e-01
124
+5.000000000000000000e-01
125
+0.000000000000000000e+00
126
+0.000000000000000000e+00
127
+0.000000000000000000e+00
128
+0.000000000000000000e+00
129
+0.000000000000000000e+00
130
+5.000000000000000000e-01
131
+5.000000000000000000e-01
132
+5.000000000000000000e-01
133
+5.000000000000000000e-01
134
+5.000000000000000000e-01
135
+0.000000000000000000e+00
136
+0.000000000000000000e+00
137
+0.000000000000000000e+00
138
+0.000000000000000000e+00
139
+0.000000000000000000e+00
140
+5.000000000000000000e-01
141
+5.000000000000000000e-01
142
+5.000000000000000000e-01
143
+5.000000000000000000e-01
144
+5.000000000000000000e-01
145
+5.000000000000000000e-01
146
+0.000000000000000000e+00
147
+0.000000000000000000e+00
148
+0.000000000000000000e+00
149
+0.000000000000000000e+00
150
+0.000000000000000000e+00
151
+5.000000000000000000e-01
152
+5.000000000000000000e-01
153
+5.000000000000000000e-01
154
+5.000000000000000000e-01
155
+5.000000000000000000e-01
156
+0.000000000000000000e+00
157
+0.000000000000000000e+00
158
+0.000000000000000000e+00
159
+0.000000000000000000e+00
160
+0.000000000000000000e+00
161
+5.000000000000000000e-01
162
+5.000000000000000000e-01
163
+5.000000000000000000e-01
164
+5.000000000000000000e-01
165
+5.000000000000000000e-01
166
+0.000000000000000000e+00
167
+0.000000000000000000e+00
168
+0.000000000000000000e+00
169
+0.000000000000000000e+00
170
+0.000000000000000000e+00
171
+5.000000000000000000e-01
172
+5.000000000000000000e-01
173
+5.000000000000000000e-01
174
+5.000000000000000000e-01
175
+5.000000000000000000e-01
176
+0.000000000000000000e+00
177
+0.000000000000000000e+00
178
+0.000000000000000000e+00
179
+0.000000000000000000e+00
180
+0.000000000000000000e+00
181
+5.000000000000000000e-01
182
+5.000000000000000000e-01
183
+5.000000000000000000e-01
184
+5.000000000000000000e-01
185
+5.000000000000000000e-01
186
+0.000000000000000000e+00
187
+0.000000000000000000e+00
188
+0.000000000000000000e+00
189
+0.000000000000000000e+00
190
+0.000000000000000000e+00
191
+5.000000000000000000e-01
192
+5.000000000000000000e-01
193
+5.000000000000000000e-01
194
+5.000000000000000000e-01
195
+5.000000000000000000e-01
196
+0.000000000000000000e+00
197
+0.000000000000000000e+00
198
+0.000000000000000000e+00
199
+0.000000000000000000e+00
200
+0.000000000000000000e+00
201
+5.000000000000000000e-01
202
+5.000000000000000000e-01

+ 102
- 0
csv/L.csv View File

@@ -0,0 +1,102 @@
1
+0.000000000000000000e+00,1.000000000000000000e+00
2
+2.000000000000000111e-01,1.000000000000000000e+00
3
+4.000000000000000222e-01,1.000000000000000000e+00
4
+6.000000000000000888e-01,1.000000000000000000e+00
5
+8.000000000000000444e-01,3.000000000000000000e+00
6
+1.000000000000000000e+00,2.000000000000000000e+00
7
+1.200000000000000178e+00,2.000000000000000000e+00
8
+1.400000000000000133e+00,2.000000000000000000e+00
9
+1.600000000000000089e+00,2.000000000000000000e+00
10
+1.800000000000000044e+00,3.000000000000000000e+00
11
+2.000000000000000000e+00,1.000000000000000000e+00
12
+2.200000000000000178e+00,1.000000000000000000e+00
13
+2.400000000000000355e+00,1.000000000000000000e+00
14
+2.600000000000000089e+00,1.000000000000000000e+00
15
+2.800000000000000266e+00,2.000000000000000000e+00
16
+3.000000000000000000e+00,2.000000000000000000e+00
17
+3.200000000000000178e+00,2.000000000000000000e+00
18
+3.400000000000000355e+00,2.000000000000000000e+00
19
+3.600000000000000089e+00,2.000000000000000000e+00
20
+3.800000000000000266e+00,1.000000000000000000e+00
21
+4.000000000000000000e+00,1.000000000000000000e+00
22
+4.200000000000000178e+00,1.000000000000000000e+00
23
+4.400000000000000355e+00,1.000000000000000000e+00
24
+4.600000000000000533e+00,1.000000000000000000e+00
25
+4.800000000000000711e+00,2.000000000000000000e+00
26
+5.000000000000000000e+00,2.000000000000000000e+00
27
+5.200000000000000178e+00,2.000000000000000000e+00
28
+5.400000000000000355e+00,2.000000000000000000e+00
29
+5.600000000000000533e+00,2.000000000000000000e+00
30
+5.800000000000000711e+00,1.000000000000000000e+00
31
+6.000000000000000000e+00,1.000000000000000000e+00
32
+6.200000000000000178e+00,1.000000000000000000e+00
33
+6.400000000000000355e+00,1.000000000000000000e+00
34
+6.600000000000000533e+00,1.000000000000000000e+00
35
+6.800000000000000711e+00,2.000000000000000000e+00
36
+7.000000000000000000e+00,2.000000000000000000e+00
37
+7.200000000000000178e+00,2.000000000000000000e+00
38
+7.400000000000000355e+00,2.000000000000000000e+00
39
+7.600000000000000533e+00,2.000000000000000000e+00
40
+7.800000000000000711e+00,1.000000000000000000e+00
41
+8.000000000000000000e+00,1.000000000000000000e+00
42
+8.200000000000001066e+00,1.000000000000000000e+00
43
+8.400000000000000355e+00,1.000000000000000000e+00
44
+8.599999999999999645e+00,1.000000000000000000e+00
45
+8.800000000000000711e+00,3.000000000000000000e+00
46
+9.000000000000000000e+00,2.000000000000000000e+00
47
+9.200000000000001066e+00,2.000000000000000000e+00
48
+9.400000000000000355e+00,2.000000000000000000e+00
49
+9.600000000000001421e+00,2.000000000000000000e+00
50
+9.800000000000000711e+00,3.000000000000000000e+00
51
+1.000000000000000000e+01,1.000000000000000000e+00
52
+1.020000000000000107e+01,1.000000000000000000e+00
53
+1.040000000000000036e+01,1.000000000000000000e+00
54
+1.060000000000000142e+01,1.000000000000000000e+00
55
+1.080000000000000071e+01,3.000000000000000000e+00
56
+1.100000000000000000e+01,2.000000000000000000e+00
57
+1.120000000000000107e+01,2.000000000000000000e+00
58
+1.140000000000000036e+01,2.000000000000000000e+00
59
+1.160000000000000142e+01,2.000000000000000000e+00
60
+1.180000000000000071e+01,3.000000000000000000e+00
61
+1.200000000000000000e+01,1.000000000000000000e+00
62
+1.220000000000000107e+01,1.000000000000000000e+00
63
+1.240000000000000036e+01,1.000000000000000000e+00
64
+1.260000000000000142e+01,1.000000000000000000e+00
65
+1.280000000000000071e+01,3.000000000000000000e+00
66
+1.300000000000000000e+01,2.000000000000000000e+00
67
+1.320000000000000107e+01,2.000000000000000000e+00
68
+1.340000000000000036e+01,2.000000000000000000e+00
69
+1.360000000000000142e+01,2.000000000000000000e+00
70
+1.380000000000000071e+01,3.000000000000000000e+00
71
+1.400000000000000000e+01,1.000000000000000000e+00
72
+1.420000000000000107e+01,1.000000000000000000e+00
73
+1.440000000000000036e+01,1.000000000000000000e+00
74
+1.460000000000000142e+01,1.000000000000000000e+00
75
+1.480000000000000071e+01,3.000000000000000000e+00
76
+1.500000000000000000e+01,2.000000000000000000e+00
77
+1.520000000000000107e+01,2.000000000000000000e+00
78
+1.540000000000000036e+01,2.000000000000000000e+00
79
+1.560000000000000142e+01,2.000000000000000000e+00
80
+1.580000000000000071e+01,3.000000000000000000e+00
81
+1.600000000000000000e+01,1.000000000000000000e+00
82
+1.619999999999999929e+01,1.000000000000000000e+00
83
+1.640000000000000213e+01,1.000000000000000000e+00
84
+1.660000000000000142e+01,1.000000000000000000e+00
85
+1.680000000000000071e+01,3.000000000000000000e+00
86
+1.700000000000000000e+01,2.000000000000000000e+00
87
+1.719999999999999929e+01,2.000000000000000000e+00
88
+1.740000000000000213e+01,2.000000000000000000e+00
89
+1.760000000000000142e+01,2.000000000000000000e+00
90
+1.780000000000000071e+01,3.000000000000000000e+00
91
+1.800000000000000000e+01,1.000000000000000000e+00
92
+1.819999999999999929e+01,1.000000000000000000e+00
93
+1.840000000000000213e+01,1.000000000000000000e+00
94
+1.860000000000000142e+01,1.000000000000000000e+00
95
+1.880000000000000071e+01,3.000000000000000000e+00
96
+1.900000000000000000e+01,2.000000000000000000e+00
97
+1.920000000000000284e+01,2.000000000000000000e+00
98
+1.940000000000000213e+01,2.000000000000000000e+00
99
+1.960000000000000142e+01,2.000000000000000000e+00
100
+1.980000000000000071e+01,3.000000000000000000e+00
101
+1.989999999999999858e+01,1.000000000000000000e+00
102
+2.000000000000000000e+01,3.000000000000000000e+00

+ 101
- 0
csv/SCAL.csv View File

@@ -0,0 +1,101 @@
1
+9.875527027155333801e-03
2
+9.880415406702709735e-03
3
+9.885205023681455441e-03
4
+9.889895873937569123e-03
5
+9.894487962810083817e-03
6
+9.898981286120364653e-03
7
+9.903375849182693633e-03
8
+9.907671647793802813e-03
9
+9.911868687243228365e-03
10
+9.915966963303065804e-03
11
+9.919966481238108938e-03
12
+9.923867236795811536e-03
13
+9.927669235216230251e-03
14
+9.931372472222175368e-03
15
+9.934976953028966384e-03
16
+9.938482673334771839e-03
17
+9.941889638330179274e-03
18
+9.945197843688710279e-03
19
+9.948407294576213505e-03
20
+9.951517986641567060e-03
21
+9.954529925025892845e-03
22
+9.957443105353416812e-03
23
+9.960257532740535852e-03
24
+9.962973202786820290e-03
25
+9.965590120583945472e-03
26
+9.968108281706824364e-03
27
+9.970527691222415972e-03
28
+9.972848344680970700e-03
29
+9.975070247124731213e-03
30
+9.977193394079288818e-03
31
+9.979217790562169843e-03
32
+9.981143432074301297e-03
33
+9.982970323608501839e-03
34
+9.984698460641034182e-03
35
+9.986327848139995439e-03
36
+9.987858481556996168e-03
37
+9.989290365835417143e-03
38
+9.990623496402199419e-03
39
+9.991857878176019572e-03
40
+9.992993506559149156e-03
41
+9.994030386445560549e-03
42
+9.994968513212850597e-03
43
+9.995807891730292688e-03
44
+9.996548517350808960e-03
45
+9.997190394918975542e-03
46
+9.997733519763025989e-03
47
+9.998177896702854783e-03
48
+9.998523521042005160e-03
49
+9.998770397575677815e-03
50
+9.998918521582748214e-03
51
+9.998967897833700713e-03
52
+9.998918521582748214e-03
53
+9.998770397575674346e-03
54
+9.998523521042012099e-03
55
+9.998177896702849579e-03
56
+9.997733519763029458e-03
57
+9.997190394918975542e-03
58
+9.996548517350808960e-03
59
+9.995807891730296157e-03
60
+9.994968513212847128e-03
61
+9.994030386445560549e-03
62
+9.992993506559142217e-03
63
+9.991857878176017838e-03
64
+9.990623496402192480e-03
65
+9.989290365835417143e-03
66
+9.987858481556992699e-03
67
+9.986327848139997174e-03
68
+9.984698460641030712e-03
69
+9.982970323608501839e-03
70
+9.981143432074301297e-03
71
+9.979217790562164639e-03
72
+9.977193394079283614e-03
73
+9.975070247124727743e-03
74
+9.972848344680962027e-03
75
+9.970527691222415972e-03
76
+9.968108281706815690e-03
77
+9.965590120583945472e-03
78
+9.962973202786811616e-03
79
+9.960257532740535852e-03
80
+9.957443105353406404e-03
81
+9.954529925025892845e-03
82
+9.951517986641556651e-03
83
+9.948407294576208301e-03
84
+9.945197843688701606e-03
85
+9.941889638330175805e-03
86
+9.938482673334764900e-03
87
+9.934976953028966384e-03
88
+9.931372472222164960e-03
89
+9.927669235216226781e-03
90
+9.923867236795802863e-03
91
+9.919966481238107203e-03
92
+9.915966963303053661e-03
93
+9.911868687243231835e-03
94
+9.907671647793790670e-03
95
+9.903375849182700572e-03
96
+9.898981286120350775e-03
97
+9.894487962810090756e-03
98
+9.889895873937553511e-03
99
+9.885205023681464115e-03
100
+9.880415406702690653e-03
101
+9.875527027155345944e-03

+ 101
- 0
csv/SF.csv View File

@@ -0,0 +1,101 @@
1
+8.279162354307141314e-01
2
+9.058144003547491385e-01
3
+9.814420386864097789e-01
4
+1.044457443180543565e+00
5
+1.108430775361605392e+00
6
+1.164536576040450244e+00
7
+1.229874438537713166e+00
8
+1.281545020903052556e+00
9
+1.337296461900192579e+00
10
+1.379213218077901981e+00
11
+1.427232463527305484e+00
12
+1.459554218814571813e+00
13
+1.494493903590228046e+00
14
+1.531075786431060237e+00
15
+1.556798002169393014e+00
16
+1.591352936621466041e+00
17
+1.615995398037212594e+00
18
+1.644922185135466508e+00
19
+1.668604691365689119e+00
20
+1.690979941635617934e+00
21
+1.713207623984341677e+00
22
+1.729528718755995342e+00
23
+1.747946242466788647e+00
24
+1.762520152861067002e+00
25
+1.778275902164144995e+00
26
+1.791827891671982087e+00
27
+1.805400699794008856e+00
28
+1.818005775943044799e+00
29
+1.828990549367728091e+00
30
+1.840217858003180051e+00
31
+1.850505289593356029e+00
32
+1.858488751029768959e+00
33
+1.866857268445507589e+00
34
+1.873739428495174053e+00
35
+1.880834871738611636e+00
36
+1.886876321823361069e+00
37
+1.892805461201894746e+00
38
+1.898605159383951246e+00
39
+1.903342947086712034e+00
40
+1.908435147992229020e+00
41
+1.912887663762121715e+00
42
+1.915572644024085447e+00
43
+1.918827909617829786e+00
44
+1.920587244825601614e+00
45
+1.923326515786471402e+00
46
+1.924179858516978925e+00
47
+1.926060902995682556e+00
48
+1.927906162135200718e+00
49
+1.928397901714244300e+00
50
+1.929888687777889800e+00
51
+1.929501187882246116e+00
52
+1.929581336838748218e+00
53
+1.928203273345603819e+00
54
+1.926777725272750486e+00
55
+1.925252666860237882e+00
56
+1.922405238411911910e+00
57
+1.920808625310032713e+00
58
+1.918757825977825959e+00
59
+1.916027340556965441e+00
60
+1.914371296486949259e+00
61
+1.909876596817369432e+00
62
+1.907154985037915917e+00
63
+1.900827044667468879e+00
64
+1.895761684486172882e+00
65
+1.889203455427365785e+00
66
+1.881222296836722041e+00
67
+1.875338851804557505e+00
68
+1.868347199244248191e+00
69
+1.861017069270422652e+00
70
+1.855355020282099421e+00
71
+1.844487976503920379e+00
72
+1.837560354859859801e+00
73
+1.823086749031125553e+00
74
+1.811798284625304900e+00
75
+1.796870691741520476e+00
76
+1.779720438165997454e+00
77
+1.766560242167598771e+00
78
+1.751073652280282555e+00
79
+1.734868994957039945e+00
80
+1.722737896145584147e+00
81
+1.699588700333208191e+00
82
+1.685184329229693079e+00
83
+1.658710560236359655e+00
84
+1.636458690026223461e+00
85
+1.605379025347835498e+00
86
+1.578060114465184682e+00
87
+1.544194887648407954e+00
88
+1.517228981872852156e+00
89
+1.488168035957224911e+00
90
+1.451838810683560244e+00
91
+1.418525012287912501e+00
92
+1.370527443091987596e+00
93
+1.328611040635534923e+00
94
+1.269800207621274435e+00
95
+1.211472571374244467e+00
96
+1.151287784965064009e+00
97
+1.086001586877506586e+00
98
+1.035757101041327788e+00
99
+9.704514977999247893e-01
100
+9.025006727015659758e-01
101
+8.274087973889836523e-01

+ 200
- 0
csv/X.csv View File

@@ -0,0 +1,200 @@
1
+8.237236283639602918e-01
2
+8.472821726185544300e-01
3
+8.647913047535317199e-01
4
+8.821818349372271273e-01
5
+8.950422855233854946e-01
6
+8.950422855233857167e-01
7
+8.950422855233857167e-01
8
+8.950422855233858277e-01
9
+8.950422855233857167e-01
10
+8.950422855233858277e-01
11
+9.084492366422537435e-01
12
+9.231958808852069875e-01
13
+9.265993846872848660e-01
14
+9.452051467178449462e-01
15
+9.388774835072988267e-01
16
+9.452051467178449462e-01
17
+9.388774835072988267e-01
18
+9.452051467178449462e-01
19
+9.388774835072988267e-01
20
+9.452051467178452793e-01
21
+9.511555823273126764e-01
22
+9.584472360253387668e-01
23
+9.595261277335117578e-01
24
+9.683651042148593291e-01
25
+9.646587633583918464e-01
26
+9.683651042148598842e-01
27
+9.646587633583917354e-01
28
+9.683651042148598842e-01
29
+9.646587633583918464e-01
30
+9.683651042148597732e-01
31
+9.697913989832721571e-01
32
+9.732213890787574373e-01
33
+9.725073465268140316e-01
34
+9.759353850689560517e-01
35
+9.730238817925007666e-01
36
+9.759353850689561627e-01
37
+9.730238817925004335e-01
38
+9.759353850689561627e-01
39
+9.730238817925003225e-01
40
+9.759353850689561627e-01
41
+9.735404170581873906e-01
42
+9.747990426701877986e-01
43
+9.719401856885295787e-01
44
+9.716466236850345162e-01
45
+9.680951691739537734e-01
46
+9.680951691739538845e-01
47
+9.680951691739539955e-01
48
+9.680951691739538845e-01
49
+9.680951691739539955e-01
50
+9.680951691739538845e-01
51
+9.650510124632859510e-01
52
+9.613392314216322454e-01
53
+9.568344043699981150e-01
54
+9.514904321830410794e-01
55
+9.451645486263099416e-01
56
+9.451645486263099416e-01
57
+9.451645486263100526e-01
58
+9.451645486263100526e-01
59
+9.451645486263100526e-01
60
+9.451645486263100526e-01
61
+9.377144708205034096e-01
62
+9.290938020670312936e-01
63
+9.190020638221960692e-01
64
+9.073505596731151757e-01
65
+8.938098219296640723e-01
66
+8.938098219296639613e-01
67
+8.938098219296639613e-01
68
+8.938098219296639613e-01
69
+8.938098219296640723e-01
70
+8.938098219296641833e-01
71
+8.780920548282996618e-01
72
+8.600939511442725616e-01
73
+8.391842215944969041e-01
74
+8.151855964504228780e-01
75
+7.874111260882538099e-01
76
+7.874111260882539209e-01
77
+7.874111260882540320e-01
78
+7.874111260882540320e-01
79
+7.874111260882540320e-01
80
+7.874111260882540320e-01
81
+7.552783221392910251e-01
82
+7.185734172702177291e-01
83
+6.760070809326180763e-01
84
+6.272215818337982807e-01
85
+5.708164062005556261e-01
86
+5.708164062005557371e-01
87
+5.708164062005557371e-01
88
+5.708164062005557371e-01
89
+5.708164062005557371e-01
90
+5.708164062005557371e-01
91
+5.056124386204988719e-01
92
+4.311746447700428586e-01
93
+3.448873248204227249e-01
94
+2.460268549211316047e-01
95
+1.317531970059818425e-01
96
+1.317531970059818702e-01
97
+1.317531970059818702e-01
98
+1.317531970059818702e-01
99
+1.317531970059818702e-01
100
+1.317531970059818702e-01
101
+1.311697539599515072e-01
102
+2.453772910738699631e-01
103
+3.440356971131438013e-01
104
+4.303847654336496187e-01
105
+4.303847654336496187e-01
106
+4.303847654336497297e-01
107
+4.303847654336496742e-01
108
+4.303847654336497297e-01
109
+4.303847654336496742e-01
110
+4.438117563085469119e-01
111
+5.880861622053110693e-01
112
+5.597054957354743454e-01
113
+7.128344519533966750e-01
114
+6.403756748212394712e-01
115
+7.128344519533961199e-01
116
+6.403756748212392491e-01
117
+7.128344519533961199e-01
118
+6.403756748212394712e-01
119
+7.128344519533965640e-01
120
+7.210458539070050410e-01
121
+7.916362540140261528e-01
122
+7.793997013053314271e-01
123
+8.537689563957779892e-01
124
+8.201055248080839633e-01
125
+8.537689563957778782e-01
126
+8.201055248080834081e-01
127
+8.537689563957778782e-01
128
+8.201055248080834081e-01
129
+8.537689563957775452e-01
130
+8.608113483108363884e-01
131
+8.925046918008538999e-01
132
+8.903820796784561731e-01
133
+9.226408025499989174e-01
134
+9.111833774203520475e-01
135
+9.226408025499990284e-01
136
+9.111833774203517144e-01
137
+9.226408025499990284e-01
138
+9.111833774203517144e-01
139
+9.226408025499991394e-01
140
+9.319846751622474779e-01
141
+9.403994417071440681e-01
142
+9.473447469171228930e-01
143
+9.533900362008603802e-01
144
+9.584923984253679308e-01
145
+9.584923984253679308e-01
146
+9.584923984253679308e-01
147
+9.584923984253679308e-01
148
+9.584923984253679308e-01
149
+9.584923984253679308e-01
150
+9.627124149130731334e-01
151
+9.661837751977747635e-01
152
+9.689912054959040066e-01
153
+9.711698539860037505e-01
154
+9.727892925184066231e-01
155
+9.727892925184066231e-01
156
+9.727892925184065120e-01
157
+9.727892925184066231e-01
158
+9.727892925184066231e-01
159
+9.727892925184066231e-01
160
+9.738508780542949328e-01
161
+9.744054761841683954e-01
162
+9.744386154057476102e-01
163
+9.739740979446640523e-01
164
+9.729814419896590794e-01
165
+9.729814419896594124e-01
166
+9.729814419896591904e-01
167
+9.729814419896594124e-01
168
+9.729814419896591904e-01
169
+9.729814419896593014e-01
170
+9.714348723388127604e-01
171
+9.693601857708454839e-01
172
+9.666362117219863270e-01
173
+9.632877444136985767e-01
174
+9.591684480429191195e-01
175
+9.591684480429193416e-01
176
+9.591684480429193416e-01
177
+9.591684480429193416e-01
178
+9.591684480429193416e-01
179
+9.591684480429194526e-01
180
+9.542120545982384350e-01
181
+9.484326638773770135e-01
182
+9.415395971620468263e-01
183
+9.335714928220242292e-01
184
+9.241903074988199185e-01
185
+9.241903074988198075e-01
186
+9.241903074988200295e-01
187
+9.241903074988198075e-01
188
+9.241903074988198075e-01
189
+9.241903074988198075e-01
190
+9.132549301666330610e-01
191
+9.007750723108800539e-01
192
+8.861398729503729976e-01
193
+8.694218429078111043e-01
194
+8.499160055701423522e-01
195
+8.499160055701422412e-01
196
+8.499160055701422412e-01
197
+8.499160055701422412e-01
198
+8.499160055701422412e-01
199
+8.499160055701422412e-01
200
+8.274037375231615421e-01

+ 202
- 0
csv/Y.csv View File

@@ -0,0 +1,202 @@
1
+8.237236283639602918e-01
2
+8.472821726185544300e-01
3
+8.647913047535317199e-01
4
+8.821818349372271273e-01
5
+8.950422855233854946e-01
6
+8.950422855233857167e-01
7
+8.950422855233857167e-01
8
+8.950422855233858277e-01
9
+8.950422855233857167e-01
10
+8.950422855233858277e-01
11
+9.084492366422537435e-01
12
+9.231958808852069875e-01
13
+9.265993846872848660e-01
14
+9.452051467178449462e-01
15
+9.388774835072988267e-01
16
+9.452051467178449462e-01
17
+9.388774835072988267e-01
18
+9.452051467178449462e-01
19
+9.388774835072988267e-01
20
+9.452051467178452793e-01
21
+9.511555823273126764e-01
22
+9.584472360253387668e-01
23
+9.595261277335117578e-01
24
+9.683651042148593291e-01
25
+9.646587633583918464e-01
26
+9.683651042148598842e-01
27
+9.646587633583917354e-01
28
+9.683651042148598842e-01
29
+9.646587633583918464e-01
30
+9.683651042148597732e-01
31
+9.697913989832721571e-01
32
+9.732213890787574373e-01
33
+9.725073465268140316e-01
34
+9.759353850689560517e-01
35
+9.730238817925007666e-01
36
+9.759353850689561627e-01
37
+9.730238817925004335e-01
38
+9.759353850689561627e-01
39
+9.730238817925003225e-01
40
+9.759353850689561627e-01
41
+9.735404170581873906e-01
42
+9.747990426701877986e-01
43
+9.719401856885295787e-01
44
+9.716466236850345162e-01
45
+9.680951691739537734e-01
46
+9.680951691739538845e-01
47
+9.680951691739539955e-01
48
+9.680951691739538845e-01
49
+9.680951691739539955e-01
50
+9.680951691739538845e-01
51
+9.650510124632859510e-01
52
+9.613392314216322454e-01
53
+9.568344043699981150e-01
54
+9.514904321830410794e-01
55
+9.451645486263099416e-01
56
+9.451645486263099416e-01
57
+9.451645486263100526e-01
58
+9.451645486263100526e-01
59
+9.451645486263100526e-01
60
+9.451645486263100526e-01
61
+9.377144708205034096e-01
62
+9.290938020670312936e-01
63
+9.190020638221960692e-01
64
+9.073505596731151757e-01
65
+8.938098219296640723e-01
66
+8.938098219296639613e-01
67
+8.938098219296639613e-01
68
+8.938098219296639613e-01
69
+8.938098219296640723e-01
70
+8.938098219296641833e-01
71
+8.780920548282996618e-01
72
+8.600939511442725616e-01
73
+8.391842215944969041e-01
74
+8.151855964504228780e-01
75
+7.874111260882538099e-01
76
+7.874111260882539209e-01
77
+7.874111260882540320e-01
78
+7.874111260882540320e-01
79
+7.874111260882540320e-01
80
+7.874111260882540320e-01
81
+7.552783221392910251e-01
82
+7.185734172702177291e-01
83
+6.760070809326180763e-01
84
+6.272215818337982807e-01
85
+5.708164062005556261e-01
86
+5.708164062005557371e-01
87
+5.708164062005557371e-01
88
+5.708164062005557371e-01
89
+5.708164062005557371e-01
90
+5.708164062005557371e-01
91
+5.056124386204988719e-01
92
+4.311746447700428586e-01
93
+3.448873248204227249e-01
94
+2.460268549211316047e-01
95
+1.317531970059818425e-01
96
+1.317531970059818702e-01
97
+1.317531970059818702e-01
98
+1.317531970059818702e-01
99
+1.317531970059818702e-01
100
+1.317531970059818702e-01
101
+0.000000000000000000e+00
102
+0.000000000000000000e+00
103
+1.311697539599515072e-01
104
+2.453772910738699631e-01
105
+3.440356971131438013e-01
106
+4.303847654336496187e-01
107
+4.303847654336496187e-01
108
+4.303847654336497297e-01
109
+4.303847654336496742e-01
110
+4.303847654336497297e-01
111
+4.303847654336496742e-01
112
+4.438117563085469119e-01
113
+5.880861622053110693e-01
114
+5.597054957354743454e-01
115
+7.128344519533966750e-01
116
+6.403756748212394712e-01
117
+7.128344519533961199e-01
118
+6.403756748212392491e-01
119
+7.128344519533961199e-01
120
+6.403756748212394712e-01
121
+7.128344519533965640e-01
122
+7.210458539070050410e-01
123
+7.916362540140261528e-01
124
+7.793997013053314271e-01
125
+8.537689563957779892e-01
126
+8.201055248080839633e-01
127
+8.537689563957778782e-01
128
+8.201055248080834081e-01
129
+8.537689563957778782e-01
130
+8.201055248080834081e-01
131
+8.537689563957775452e-01
132
+8.608113483108363884e-01
133
+8.925046918008538999e-01
134
+8.903820796784561731e-01
135
+9.226408025499989174e-01
136
+9.111833774203520475e-01
137
+9.226408025499990284e-01
138
+9.111833774203517144e-01
139
+9.226408025499990284e-01
140
+9.111833774203517144e-01
141
+9.226408025499991394e-01
142
+9.319846751622474779e-01
143
+9.403994417071440681e-01
144
+9.473447469171228930e-01
145
+9.533900362008603802e-01
146
+9.584923984253679308e-01
147
+9.584923984253679308e-01
148
+9.584923984253679308e-01
149
+9.584923984253679308e-01
150
+9.584923984253679308e-01
151
+9.584923984253679308e-01
152
+9.627124149130731334e-01
153
+9.661837751977747635e-01
154
+9.689912054959040066e-01
155
+9.711698539860037505e-01
156
+9.727892925184066231e-01
157
+9.727892925184066231e-01
158
+9.727892925184065120e-01
159
+9.727892925184066231e-01
160
+9.727892925184066231e-01
161
+9.727892925184066231e-01
162
+9.738508780542949328e-01
163
+9.744054761841683954e-01
164
+9.744386154057476102e-01
165
+9.739740979446640523e-01
166
+9.729814419896590794e-01
167
+9.729814419896594124e-01
168
+9.729814419896591904e-01
169
+9.729814419896594124e-01
170
+9.729814419896591904e-01
171
+9.729814419896593014e-01
172
+9.714348723388127604e-01
173
+9.693601857708454839e-01
174
+9.666362117219863270e-01
175
+9.632877444136985767e-01
176
+9.591684480429191195e-01
177
+9.591684480429193416e-01
178
+9.591684480429193416e-01
179
+9.591684480429193416e-01
180
+9.591684480429193416e-01
181
+9.591684480429194526e-01
182
+9.542120545982384350e-01
183
+9.484326638773770135e-01
184
+9.415395971620468263e-01
185
+9.335714928220242292e-01
186
+9.241903074988199185e-01
187
+9.241903074988198075e-01
188
+9.241903074988200295e-01
189
+9.241903074988198075e-01
190
+9.241903074988198075e-01
191
+9.241903074988198075e-01
192
+9.132549301666330610e-01
193
+9.007750723108800539e-01
194
+8.861398729503729976e-01
195
+8.694218429078111043e-01
196
+8.499160055701423522e-01
197
+8.499160055701422412e-01
198
+8.499160055701422412e-01
199
+8.499160055701422412e-01
200
+8.499160055701422412e-01
201
+8.499160055701422412e-01
202
+8.274037375231615421e-01

+ 202
- 0
csv/Z.csv View File

@@ -0,0 +1,202 @@
1
+8.237236283639602918e-01
2
+8.472821726185544300e-01
3
+8.647913047535317199e-01
4
+8.821818349372271273e-01
5
+8.950422855233854946e-01
6
+8.950422855233857167e-01
7
+8.950422855233857167e-01
8
+8.950422855233858277e-01
9
+8.950422855233857167e-01
10
+8.950422855233858277e-01
11
+9.084492366422537435e-01
12
+9.231958808852069875e-01
13
+9.265993846872848660e-01
14
+9.452051467178449462e-01
15
+9.388774835072988267e-01
16
+9.452051467178449462e-01
17
+9.388774835072988267e-01
18
+9.452051467178449462e-01
19
+9.388774835072988267e-01
20
+9.452051467178452793e-01
21
+9.511555823273126764e-01
22
+9.584472360253387668e-01
23
+9.595261277335117578e-01
24
+9.683651042148593291e-01
25
+9.646587633583918464e-01
26
+9.683651042148598842e-01
27
+9.646587633583917354e-01
28
+9.683651042148598842e-01
29
+9.646587633583918464e-01
30
+9.683651042148597732e-01
31
+9.697913989832721571e-01
32
+9.732213890787574373e-01
33
+9.725073465268140316e-01
34
+9.759353850689560517e-01
35
+9.730238817925007666e-01
36
+9.759353850689561627e-01
37
+9.730238817925004335e-01
38
+9.759353850689561627e-01
39
+9.730238817925003225e-01
40
+9.759353850689561627e-01
41
+9.735404170581873906e-01
42
+9.747990426701877986e-01
43
+9.719401856885295787e-01
44
+9.716466236850345162e-01
45
+9.680951691739537734e-01
46
+9.680951691739538845e-01
47
+9.680951691739539955e-01
48
+9.680951691739538845e-01
49
+9.680951691739539955e-01
50
+9.680951691739538845e-01
51
+9.650510124632859510e-01
52
+9.613392314216322454e-01
53
+9.568344043699981150e-01
54
+9.514904321830410794e-01
55
+9.451645486263099416e-01
56
+9.451645486263099416e-01
57
+9.451645486263100526e-01
58
+9.451645486263100526e-01
59
+9.451645486263100526e-01
60
+9.451645486263100526e-01
61
+9.377144708205034096e-01
62
+9.290938020670312936e-01
63
+9.190020638221960692e-01
64
+9.073505596731151757e-01
65
+8.938098219296640723e-01
66
+8.938098219296639613e-01
67
+8.938098219296639613e-01
68
+8.938098219296639613e-01
69
+8.938098219296640723e-01
70
+8.938098219296641833e-01
71
+8.780920548282996618e-01
72
+8.600939511442725616e-01
73
+8.391842215944969041e-01
74
+8.151855964504228780e-01
75
+7.874111260882538099e-01
76
+7.874111260882539209e-01
77
+7.874111260882540320e-01
78
+7.874111260882540320e-01
79
+7.874111260882540320e-01
80
+7.874111260882540320e-01
81
+7.552783221392910251e-01
82
+7.185734172702177291e-01
83
+6.760070809326180763e-01
84
+6.272215818337982807e-01
85
+5.708164062005556261e-01
86
+5.708164062005557371e-01
87
+5.708164062005557371e-01
88
+5.708164062005557371e-01
89
+5.708164062005557371e-01
90
+5.708164062005557371e-01
91
+5.056124386204988719e-01
92
+4.311746447700428586e-01
93
+3.448873248204227249e-01
94
+2.460268549211316047e-01
95
+1.317531970059818425e-01
96
+1.317531970059818702e-01
97
+1.317531970059818702e-01
98
+1.317531970059818702e-01
99
+1.317531970059818702e-01
100
+1.317531970059818702e-01
101
+6.813169257258520906e-02
102
+1.311697539599515072e-01
103
+2.453772910738699631e-01
104
+3.440356971131438013e-01
105
+4.303847654336496187e-01
106
+4.303847654336496187e-01
107
+4.303847654336497297e-01
108
+4.303847654336496742e-01
109
+4.303847654336497297e-01
110
+4.303847654336496742e-01
111
+4.438117563085469119e-01
112
+5.880861622053110693e-01
113
+5.597054957354743454e-01
114
+7.128344519533966750e-01
115
+6.403756748212394712e-01
116
+7.128344519533961199e-01
117
+6.403756748212392491e-01
118
+7.128344519533961199e-01
119
+6.403756748212394712e-01
120
+7.128344519533965640e-01
121
+7.210458539070050410e-01
122
+7.916362540140261528e-01
123
+7.793997013053314271e-01
124
+8.537689563957779892e-01
125
+8.201055248080839633e-01
126
+8.537689563957778782e-01
127
+8.201055248080834081e-01
128
+8.537689563957778782e-01
129
+8.201055248080834081e-01
130
+8.537689563957775452e-01
131
+8.608113483108363884e-01
132
+8.925046918008538999e-01
133
+8.903820796784561731e-01
134
+9.226408025499989174e-01
135
+9.111833774203520475e-01
136
+9.226408025499990284e-01
137
+9.111833774203517144e-01
138
+9.226408025499990284e-01
139
+9.111833774203517144e-01
140
+9.226408025499991394e-01
141
+9.319846751622474779e-01
142
+9.403994417071440681e-01
143
+9.473447469171228930e-01
144
+9.533900362008603802e-01
145
+9.584923984253679308e-01
146
+9.584923984253679308e-01
147
+9.584923984253679308e-01
148
+9.584923984253679308e-01
149
+9.584923984253679308e-01
150
+9.584923984253679308e-01
151
+9.627124149130731334e-01
152
+9.661837751977747635e-01
153
+9.689912054959040066e-01
154
+9.711698539860037505e-01
155
+9.727892925184066231e-01
156
+9.727892925184066231e-01
157
+9.727892925184065120e-01
158
+9.727892925184066231e-01
159
+9.727892925184066231e-01
160
+9.727892925184066231e-01
161
+9.738508780542949328e-01
162
+9.744054761841683954e-01
163
+9.744386154057476102e-01
164
+9.739740979446640523e-01
165
+9.729814419896590794e-01
166
+9.729814419896594124e-01
167
+9.729814419896591904e-01
168
+9.729814419896594124e-01
169
+9.729814419896591904e-01
170
+9.729814419896593014e-01
171
+9.714348723388127604e-01
172
+9.693601857708454839e-01
173
+9.666362117219863270e-01
174
+9.632877444136985767e-01
175
+9.591684480429191195e-01
176
+9.591684480429193416e-01
177
+9.591684480429193416e-01
178
+9.591684480429193416e-01
179
+9.591684480429193416e-01
180
+9.591684480429194526e-01
181
+9.542120545982384350e-01
182
+9.484326638773770135e-01
183
+9.415395971620468263e-01
184
+9.335714928220242292e-01
185
+9.241903074988199185e-01
186
+9.241903074988198075e-01
187
+9.241903074988200295e-01
188
+9.241903074988198075e-01
189
+9.241903074988198075e-01
190
+9.241903074988198075e-01
191
+9.132549301666330610e-01
192
+9.007750723108800539e-01
193
+8.861398729503729976e-01
194
+8.694218429078111043e-01
195
+8.499160055701423522e-01
196
+8.499160055701422412e-01
197
+8.499160055701422412e-01
198
+8.499160055701422412e-01
199
+8.499160055701422412e-01
200
+8.499160055701422412e-01
201
+8.394590178373988465e-01
202
+8.274037375231615421e-01

+ 100
- 0
csv/extra.csv View File

@@ -0,0 +1,100 @@
1
+0.000000000000000000e+00
2
+0.000000000000000000e+00
3
+0.000000000000000000e+00
4
+0.000000000000000000e+00
5
+0.000000000000000000e+00
6
+0.000000000000000000e+00
7
+0.000000000000000000e+00
8
+0.000000000000000000e+00
9
+0.000000000000000000e+00
10
+0.000000000000000000e+00
11
+0.000000000000000000e+00
12
+0.000000000000000000e+00
13
+0.000000000000000000e+00
14
+0.000000000000000000e+00
15
+0.000000000000000000e+00
16
+0.000000000000000000e+00
17
+0.000000000000000000e+00
18
+0.000000000000000000e+00
19
+0.000000000000000000e+00
20
+0.000000000000000000e+00
21
+0.000000000000000000e+00
22
+0.000000000000000000e+00
23
+0.000000000000000000e+00
24
+0.000000000000000000e+00
25
+0.000000000000000000e+00
26
+0.000000000000000000e+00
27
+0.000000000000000000e+00
28
+0.000000000000000000e+00
29
+0.000000000000000000e+00
30
+0.000000000000000000e+00
31
+0.000000000000000000e+00
32
+0.000000000000000000e+00
33
+0.000000000000000000e+00
34
+0.000000000000000000e+00
35
+0.000000000000000000e+00
36
+0.000000000000000000e+00
37
+0.000000000000000000e+00
38
+0.000000000000000000e+00
39
+0.000000000000000000e+00
40
+0.000000000000000000e+00
41
+0.000000000000000000e+00
42
+0.000000000000000000e+00
43
+0.000000000000000000e+00
44
+0.000000000000000000e+00
45
+0.000000000000000000e+00
46
+0.000000000000000000e+00
47
+0.000000000000000000e+00
48
+0.000000000000000000e+00
49
+0.000000000000000000e+00
50
+0.000000000000000000e+00
51
+0.000000000000000000e+00
52
+0.000000000000000000e+00
53
+0.000000000000000000e+00
54
+0.000000000000000000e+00
55
+0.000000000000000000e+00
56
+0.000000000000000000e+00
57
+0.000000000000000000e+00
58
+0.000000000000000000e+00
59
+0.000000000000000000e+00
60
+0.000000000000000000e+00
61
+0.000000000000000000e+00
62
+0.000000000000000000e+00
63
+0.000000000000000000e+00
64
+0.000000000000000000e+00
65
+0.000000000000000000e+00
66
+0.000000000000000000e+00
67
+0.000000000000000000e+00
68
+0.000000000000000000e+00
69
+0.000000000000000000e+00
70
+0.000000000000000000e+00
71
+0.000000000000000000e+00
72
+0.000000000000000000e+00
73
+0.000000000000000000e+00
74
+0.000000000000000000e+00
75
+0.000000000000000000e+00
76
+0.000000000000000000e+00
77
+0.000000000000000000e+00
78
+0.000000000000000000e+00
79
+0.000000000000000000e+00
80
+0.000000000000000000e+00
81
+0.000000000000000000e+00
82
+0.000000000000000000e+00
83
+0.000000000000000000e+00
84
+0.000000000000000000e+00
85
+0.000000000000000000e+00
86
+0.000000000000000000e+00
87
+0.000000000000000000e+00
88
+0.000000000000000000e+00
89
+0.000000000000000000e+00
90
+0.000000000000000000e+00
91
+0.000000000000000000e+00
92
+0.000000000000000000e+00
93
+0.000000000000000000e+00
94
+0.000000000000000000e+00
95
+0.000000000000000000e+00
96
+0.000000000000000000e+00
97
+0.000000000000000000e+00
98
+0.000000000000000000e+00
99
+0.000000000000000000e+00
100
+1.000000000000000000e+00

+ 101
- 0
csv/h.csv View File

@@ -0,0 +1,101 @@
1
+2.000000000000000111e-01
2
+2.000000000000000111e-01
3
+2.000000000000000111e-01
4
+2.000000000000000111e-01
5
+1.999999999999999556e-01
6
+2.000000000000001776e-01
7
+1.999999999999999556e-01
8
+1.999999999999999556e-01
9
+1.999999999999999556e-01
10
+1.999999999999999556e-01
11
+2.000000000000001776e-01
12
+2.000000000000001776e-01
13
+1.999999999999997335e-01
14
+2.000000000000001776e-01
15
+1.999999999999997335e-01
16
+2.000000000000001776e-01
17
+2.000000000000001776e-01
18
+1.999999999999997335e-01
19
+2.000000000000001776e-01
20
+1.999999999999997335e-01
21
+2.000000000000001776e-01
22
+2.000000000000001776e-01
23
+2.000000000000001776e-01
24
+2.000000000000001776e-01
25
+1.999999999999992895e-01
26
+2.000000000000001776e-01
27
+2.000000000000001776e-01
28
+2.000000000000001776e-01
29
+2.000000000000001776e-01
30
+1.999999999999992895e-01
31
+2.000000000000001776e-01
32
+2.000000000000001776e-01
33
+2.000000000000001776e-01
34
+2.000000000000001776e-01
35
+1.999999999999992895e-01
36
+2.000000000000001776e-01
37
+2.000000000000001776e-01
38
+2.000000000000001776e-01
39
+2.000000000000001776e-01
40
+1.999999999999992895e-01
41
+2.000000000000010658e-01
42
+1.999999999999992895e-01
43
+1.999999999999992895e-01
44
+2.000000000000010658e-01
45
+1.999999999999992895e-01
46
+2.000000000000010658e-01
47
+1.999999999999992895e-01
48
+2.000000000000010658e-01
49
+1.999999999999992895e-01
50
+1.999999999999992895e-01
51
+2.000000000000010658e-01
52
+1.999999999999992895e-01
53
+2.000000000000010658e-01
54
+1.999999999999992895e-01
55
+1.999999999999992895e-01
56
+2.000000000000010658e-01
57
+1.999999999999992895e-01
58
+2.000000000000010658e-01
59
+1.999999999999992895e-01
60
+1.999999999999992895e-01
61
+2.000000000000010658e-01
62
+1.999999999999992895e-01
63
+2.000000000000010658e-01
64
+1.999999999999992895e-01
65
+1.999999999999992895e-01
66
+2.000000000000010658e-01
67
+1.999999999999992895e-01
68
+2.000000000000010658e-01
69
+1.999999999999992895e-01
70
+1.999999999999992895e-01
71
+2.000000000000010658e-01
72
+1.999999999999992895e-01
73
+2.000000000000010658e-01
74
+1.999999999999992895e-01
75
+1.999999999999992895e-01
76
+2.000000000000010658e-01
77
+1.999999999999992895e-01
78
+2.000000000000010658e-01
79
+1.999999999999992895e-01
80
+1.999999999999992895e-01
81
+1.999999999999992895e-01
82
+2.000000000000028422e-01
83
+1.999999999999992895e-01
84
+1.999999999999992895e-01
85
+1.999999999999992895e-01
86
+1.999999999999992895e-01
87
+2.000000000000028422e-01
88
+1.999999999999992895e-01
89
+1.999999999999992895e-01
90
+1.999999999999992895e-01
91
+1.999999999999992895e-01
92
+2.000000000000028422e-01
93
+1.999999999999992895e-01
94
+1.999999999999992895e-01
95
+1.999999999999992895e-01
96
+2.000000000000028422e-01
97
+1.999999999999992895e-01
98
+1.999999999999992895e-01
99
+1.999999999999992895e-01
100
+9.999999999999786837e-02
101
+1.000000000000014211e-01

+ 21
- 0
csv/xsn.csv View File

@@ -0,0 +1,21 @@
1
+8.000000000000000444e-01,1.000000000000000000e+00
2
+1.800000000000000044e+00,2.000000000000000000e+00
3
+2.799999999999999822e+00,1.000000000000000000e+00
4
+3.799999999999999822e+00,2.000000000000000000e+00
5
+4.799999999999999822e+00,1.000000000000000000e+00
6
+5.799999999999999822e+00,2.000000000000000000e+00
7
+6.799999999999999822e+00,1.000000000000000000e+00
8
+7.799999999999999822e+00,2.000000000000000000e+00
9
+8.800000000000000711e+00,1.000000000000000000e+00
10
+9.800000000000000711e+00,2.000000000000000000e+00
11
+1.080000000000000071e+01,1.000000000000000000e+00
12
+1.180000000000000071e+01,2.000000000000000000e+00
13
+1.280000000000000071e+01,1.000000000000000000e+00
14
+1.380000000000000071e+01,2.000000000000000000e+00
15
+1.480000000000000071e+01,1.000000000000000000e+00
16
+1.580000000000000071e+01,2.000000000000000000e+00
17
+1.680000000000000071e+01,1.000000000000000000e+00
18
+1.780000000000000071e+01,2.000000000000000000e+00
19
+1.880000000000000071e+01,1.000000000000000000e+00
20
+1.980000000000000071e+01,2.000000000000000000e+00
21
+2.000000000000000000e+01,1.000000000000000000e+00

+ 10
- 0
steady-state.config View File

@@ -0,0 +1,10 @@
1
+% Steady State Configuration file
2
+% Please use the following format for input parameters
3
+% [SLAB/ROD Geometry, # of Discrete Ordinate Directions, total length of system (value will be ignored if ROD geometry), number of points you want to calculate, total cross section, scattering ratio, boundary value in positive direction, boundary value in negative direction, homogeneous isotropic source]
4
+% ROD Example: [ROD, 0, 0.1, 100, 0.1, 0, 0, 0, 0.2]
5
+
6
+%[ROD, 0, 0.1, 100, 0.5, 0.5, 0, 0, 0.2]
7
+
8
+% periodic 
9
+% [T, n, yo, y_, m1, m2, Et1, cc, Q1, Et2, cc2, Q2]
10
+[20.0, 100, 0.0, 0.0, 1.0, 1.0, 1.0, 0.5, 1.0, 0.0, 0.0, 0.0]

+ 838
- 0
steady-state.py View File

@@ -0,0 +1,838 @@
1
+# %THIS CODE SOLVES THE STEADY-STATE, MONOENERGETIC TRANSPORT EQUATION 
2
+# %IN A HOMOGENEOUS MEDIUM WITH ISOTROPIC SCATTERING AND ISOTROPIC SOURCE
3
+# %IN ROD AND SLAB (SN) GEOMETRIES. IT USES 2-POINT CENTRAL DIFFERENCING
4
+# %(ORDER H^2) WITH 3-POINT FORWARD/BACKWARD AT THE BOUNDARIES.
5
+
6
+# %MAIN OUTPUTS: transmission "TR"; reflection "RL"
7
+# %			  scalar flux "SCAL"
8
+
9
+# %RICHARD VASQUES & NITIN KUMAR YADAV
10