PLaSK library

Loading...
Searching...
No Matches
kublybr.cc
Go to the documentation of this file.
1#include "kublybr.h"
2
3#ifndef M_PI
4#define M_PI 3.14159265358979323846
5#endif
6
7#ifdef RPSMES
8// Tu sa definicje, ktore powinny byc w C++, ale VS6 ich nie ma
9
10#include "Faddeeva.h"
11 using Faddeeva::erf;
12
13#include <stdlib.h>
14
15#define abs fabs
16
17 inline double min(double _x1, double _x2)
18{
19 if(_x1 < _x2)
20 return _x1;
21 else
22 return _x2;
23}
24
25 inline double max(double _x1, double _x2)
26{
27 if(_x1 > _x2)
28 return _x1;
29 else
30 return _x2;
31}
32
33 namespace std
34{
35 template <class IT> IT prev(IT it)
36{
37 std::advance(it, -1);
38 return it;
39}
40} // namespace std
41
42#else
43
44 using std::abs;
45 using std::max;
46 using std::min;
47
48#endif
49
50#ifdef PLASK
51#include <plask/phys/functions.hpp>
52 using plask::phys::nm_to_eV;
53#else
54 inline static double nm_to_eV(double wavelength) { return 1239.84193009 / wavelength; } // h_eV*c*1e9 = 1239.84193009
55#endif
56
57#ifndef M_PI
58#define M_PI 3.14159265358979323846
59#endif
60
61 namespace kubly
62{
63
64 const double struktura::przelm = 10 * 1.05459 / (sqrt(1.60219 * 9.10956));
65 const double struktura::przels = 1.05459 / 1.60219 * 1e-3;
66 const double struktura::kB = 1.38062 / 1.60219 * 1e-4;
67 const double struktura::eps0 = 8.8542 * 1.05459 / (100 * 1.60219 * sqrt(1.60219 * 9.10956));
68 const double struktura::c = 300 * sqrtl(9.10956 / 1.60219);
69 const double struktura::pi = M_PI;
70
71/*****************************************************************************
72 warstwa::warstwa(const warstwa & war) : x_pocz(war.x_pocz), x_kon(war.x_kon), y_pocz(war.y_pocz), y_kon(war.y_kon),
73 masa(war.masa) {}
74*****************************************************************************
75 warstwa & warstwa::operator=(const warstwa & war)
76{
77 return *this;
78}
79*****************************************************************************/
80 warstwa::warstwa(double m_p, double m_r, double x_p, double y_p, double x_k, double y_k, double niepar, double niepar2)
81 : x_pocz(x_p / struktura::przelm), x_kon(x_k / struktura::przelm), y_pocz(y_p), y_kon(y_k), nieparab(niepar),
82 nieparab_2(niepar2), m_p(m_p), nast(NULL), masa_r(m_r) // Polozenia w A
83{
84 // std::cout << "\tmasa_r w warstwa_konstruktor: " << masa_r << ", x_p: " << x_p << ", x_k: " << x_k << "\n"; // TEST
85 // // TEST LUKASZ
86 if(x_k <= x_p)
87 {
88 Error err("Zle dane!\n");
89 err << "pocz = " << x_p << "\tkoniec = " << x_k << "\tmasa_p = " << m_p << "\n";
90 throw err;
91 }
92 // std::clog<<"x_pocz = "<<x_pocz<<"\tx_kon = "<<x_kon<<"\n";
93 pole = (y_kon - y_pocz) / (x_kon - x_pocz);
94}
95/*****************************************************************************/
96 warstwa_skraj::warstwa_skraj(strona lczyp, double m_p, double m_r, double x0, double y0)
97 : warstwa(m_p, m_r, (lczyp == lewa) ? x0 - 1 : x0, y0, (lczyp == lewa) ? x0 : x0 + 1, y0), lp(lczyp), masa_p(m_p),
98 masa_r(m_r), iks(x0 / struktura::przelm), y(y0)
99{
100} // x0 w A
101/*****************************************************************************/
102 warstwa_skraj::warstwa_skraj() : warstwa(0., 0., 0., 0., 1., 0.) {}
103/*****************************************************************************/
104 warstwa_skraj::warstwa_skraj(const warstwa_skraj & war)
105 : warstwa(war.masa_p, war.masa_r, (war.lp == lewa) ? war.iks - 1 : war.iks, war.y,
106 (war.lp == lewa) ? war.iks : war.iks + 1, war.y),
107 lp(war.lp), iks(war.iks), y(war.y)
108{
109}
110/*****************************************************************************/
111//warstwa_skraj::warstwa_skraj(const warstwa_skraj & war) : warstwa(war.masa_p, war.masa_r, (war.lp == lewa)?war.iks - 1:war.iks, war.y, (war.lp == lewa)?war.iks:war.iks + 1, war.y), lp(war.lp), iks(war.iks), y(war.y) {}
112/*****************************************************************************/
113 inline double warstwa::masa_p(double E) const
114{
115 double wynik;
116 double Ek = E - (y_pocz + y_kon) / 2;
117 if((nieparab == 0 && nieparab_2 == 0) || Ek < 0)
118 {
119 wynik = m_p;
120 }
121 else
122 {
123 if((nieparab_2 < 0) && (Ek > -nieparab / (2 * nieparab_2))) // czy nie przeszlo na opadajaca czesc?
124 wynik = m_p * (1. - nieparab * nieparab / (4 * nieparab_2));
125 else
126 wynik = (1. + nieparab * Ek + nieparab_2 * Ek * Ek) * m_p;
127 }
128 return wynik;
129}
130/*****************************************************************************/
131 void warstwa::przesun_igreki(double dE)
132{
133 y_pocz += dE;
134 y_kon += dE;
135}
136/*****************************************************************************/
137 double warstwa::ffala(double x, double E) const
138{
139 double wartosc;
140 // std::clog<<" E = "<<E<<"\n";
141 if(pole != 0)
142 {
143 // std::clog<<"\n poczatek warstwy w "<<x_pocz<<" Airy";
144 // wartosc = Ai(x, E);
145 wartosc = Ai_skala(x, E);
146 }
147 else
148 {
149 if(E >= y_pocz)
150 {
151 // std::clog<<"\n poczatek warstwy w "<<x_pocz<<" tryg";
152 wartosc = tryga(x, E);
153 }
154 else
155 {
156 // std::clog<<"\n poczatek warstwy w "<<x_pocz<<" exp";
157 wartosc = expa(x, E);
158 }
159 }
160 return wartosc;
161}
162/*****************************************************************************/
163 double warstwa::ffala_prim(double x, double E) const
164{
165 double wartosc;
166 if(pole != 0)
167 {
168 // wartosc = Ai_prim(x, E);
169 wartosc = Ai_prim_skala(x, E);
170 }
171 else
172 {
173 if(E >= y_pocz)
174 {
175 wartosc = tryga_prim(x, E);
176 }
177 else
178 {
179 wartosc = expa_prim(x, E);
180 }
181 }
182 return wartosc;
183}
184/*****************************************************************************/
185 double warstwa::ffalb(double x, double E) const
186{
187 double wartosc;
188 if(pole != 0)
189 {
190 // wartosc = Bi(x, E);
191 wartosc = Bi_skala(x, E);
192 }
193 else
194 {
195 if(E >= y_pocz)
196 {
197 wartosc = trygb(x, E);
198 }
199 else
200 {
201 wartosc = expb(x, E);
202 }
203 }
204 return wartosc;
205}
206/*****************************************************************************/
207 double warstwa::ffalb_prim(double x, double E) const
208{
209 double wartosc;
210 if(pole != 0)
211 {
212 // wartosc = Bi_prim(x, E);
213 wartosc = Bi_prim_skala(x, E);
214 }
215 else
216 {
217 if(E >= y_pocz)
218 {
219 wartosc = trygb_prim(x, E);
220 }
221 else
222 {
223 wartosc = expb_prim(x, E);
224 }
225 }
226 return wartosc;
227}
228/*****************************************************************************/
229 double warstwa::tryga(double x, double E) const
230{
231 if((y_kon != y_pocz) || (E < y_pocz))
232 {
233 Error err; // MD
234 err << "Zla funkcja (tryga)!\n";
235 throw err;
236 }
237 double k = sqrt(2 * masa_p(E) * (E - y_pocz));
238 return sin(k * x);
239}
240/*****************************************************************************/
241 double warstwa::tryga_prim(double x, double E) const
242{
243 if((y_kon != y_pocz) || (E < y_pocz))
244 {
245 Error err; // MD
246 err << "Zla funkcja (tryga_prim)!\n";
247 throw err;
248 }
249 double k = sqrt(2 * masa_p(E) * (E - y_pocz));
250 return k * cos(k * x);
251}
252/*****************************************************************************/
253 double warstwa::trygb(double x, double E) const
254{
255 if((y_kon != y_pocz) || (E < y_pocz))
256 {
257 Error err; // MD
258 err << "Zla funkcja (trygb)!\n";
259 throw err;
260 }
261 double k = sqrt(2 * masa_p(E) * (E - y_pocz));
262 return cos(k * x);
263}
264/*****************************************************************************/
265 double warstwa::trygb_prim(double x, double E) const
266{
267 if((y_kon != y_pocz) || (E < y_pocz))
268 {
269 Error err; // MD
270 err << "Zla funkcja (trygb_prim)!\n";
271 throw err;
272 }
273 double k = sqrt(2 * masa_p(E) * (E - y_pocz));
274 return -k * sin(k * x);
275}
276/*****************************************************************************/
277 double warstwa::tryg_kwadr_pierwotna(double x, double E, double A, double B) const
278{
279 if((y_kon != y_pocz) || (E <= y_pocz))
280 {
281 Error err; // MD
282 err << "Zla funkcja (tryg_kwadr_pierwotna)!\n";
283 throw err;
284 }
285 double k = sqrt(2 * masa_p(E) * (E - y_pocz));
286 double si2 = sin(2 * k * x);
287 double co = cos(k * x);
288 return (A * A + B * B) * x / 2 + (si2 * (B * B - A * A) / 4 - A * B * co * co) / k;
289}
290/*****************************************************************************/
291 double warstwa::expa(double x, double E) const
292{
293 if((y_kon != y_pocz) || (E > y_pocz))
294 {
295 Error err; // MD
296 err << "Zla funkcja (expa)!\n";
297 throw err;
298 }
299 double kp = sqrt(2 * masa_p(E) * (y_pocz - E));
300 return exp(-kp * (x - x_pocz));
301}
302/*****************************************************************************/
303 double warstwa::expa_prim(double x, double E) const
304{
305 if((y_kon != y_pocz) || (E > y_pocz))
306 {
307 Error err; // MD
308 err << "Zla funkcja (expa_prim)!\n";
309 throw err;
310 }
311 double kp = sqrt(2 * masa_p(E) * (y_pocz - E));
312 return -kp * exp(-kp * (x - x_pocz));
313}
314/*****************************************************************************/
315 double warstwa::expb(double x, double E) const
316{
317 if((y_kon != y_pocz) || (E > y_pocz))
318 {
319 Error err; // MD
320 err << "Zla funkcja (expb)!\n"
321 << "y_pocz = " << y_pocz << "\ty_kon = " << y_kon << "\n";
322 throw err;
323 }
324 double kp = sqrt(2 * masa_p(E) * (y_pocz - E));
325 return exp(kp * (x - x_kon));
326}
327/*****************************************************************************/
328 double warstwa::expb_prim(double x, double E) const
329{
330 if((y_kon != y_pocz) || (E > y_pocz))
331 {
332 Error err; // MD
333 err << "Zla funkcja (expb_prim)!\n";
334 throw err;
335 }
336 double kp = sqrt(2 * masa_p(E) * (y_pocz - E));
337 return kp * exp(kp * (x - x_kon));
338}
339/*****************************************************************************/
340 double warstwa::exp_kwadr_pierwotna(double x, double E, double A, double B) const
341{
342 if((y_kon != y_pocz) || (E > y_pocz))
343 {
344 Error err; // MD
345 err << "Zla funkcja (expa)!\n";
346 throw err;
347 }
348 double kp = sqrt(2 * masa_p(E) * (y_pocz - E));
349 double b = expb(x, E);
350 double a = expa(x, E);
351 return (B * B * b * b - A * A * a * a) / (2 * kp) + 2 * A * B * x * exp(kp * (x_pocz - x_kon));
352}
353/*****************************************************************************/
354 double warstwa::Ai(double x, double E) const
355{
356 if(y_kon == y_pocz)
357 {
358 Error err; // MD
359 err << "Zla funkcja!\n";
360 throw err;
361 }
362 // rownanie: -f''(x) + (b + ax)f(x) = 0
363 // a = 2m*pole/h^2
364 // b = 2m(U - E)/h^2
365 // rozw f(x) = Ai( (ax+b)/a^{2/3} ) = Ai( a^{1/3} (x + b/a^{2/3}) )
366 double U = y_pocz - pole * x_pocz;
367 double a13 = (pole > 0) ? pow(2 * masa_p(E) * pole, 1. / 3) : -pow(-2 * masa_p(E) * pole, 1. / 3); // a^{1/3}
368 double b_a23 = (U - E) / pole;
369 double arg = a13 * (x + b_a23);
370 // std::cerr<<"\narg_Ai = "<<arg;
371 return gsl_sf_airy_Ai(arg, GSL_PREC_DOUBLE);
372}
373/*****************************************************************************/
374 double warstwa::Ai_skala(double x, double E) const
375{
376 if(y_kon == y_pocz)
377 {
378 Error err; // MD
379 err << "Zla funkcja!\n";
380 throw err;
381 }
382 // rownanie: -f''(x) + (b + ax)f(x) = 0
383 // a = 2m*pole/h^2
384 // b = 2m(U - E)/h^2
385 // rozw f(x) = Ai( (ax+b)/a^{2/3} ) = Ai( a^{1/3} (x + b/a^{2/3}) )
386 double U = y_pocz - pole * x_pocz;
387 double a13 = (pole > 0) ? pow(2 * masa_p(E) * pole, 1. / 3) : -pow(-2 * masa_p(E) * pole, 1. / 3); // a^{1/3}
388 double b_a23 = (U - E) / pole;
389 double arg = a13 * (x + b_a23);
390 double arg_pocz = a13 * (x_pocz + b_a23);
391 double przeskal = 1.0;
392 if(arg > 0 && arg_pocz > 0)
393 {
394 przeskal = exp(-2 * (pow(arg, 3. / 2) - pow(arg_pocz, 3. / 2)) / 3);
395 }
396 else
397 {
398 if(arg > 0)
399 {
400 przeskal = exp(-2 * pow(arg, 3. / 2) / 3);
401 }
402 if(arg_pocz > 0)
403 {
404 przeskal = exp(2 * pow(arg_pocz, 3. / 2) / 3);
405 }
406 }
407 // std::clog<<"arg = "<<arg<<"\targpocz = "<<arg_pocz<<"\tprzeskal = "<<przeskal<<"\n";
408 return gsl_sf_airy_Ai_scaled(arg, GSL_PREC_DOUBLE) * przeskal; // Na razie nie ma
409}
410/*****************************************************************************/
411 double warstwa::Ai_prim(double x, double E) const
412{
413 if(y_kon == y_pocz)
414 {
415 Error err; // MD
416 err << "Zla funkcja!\n";
417 throw err;
418 }
419 double U = y_pocz - pole * x_pocz;
420 double a13 = (pole > 0) ? pow(2 * masa_p(E) * pole, 1. / 3) : -pow(-2 * masa_p(E) * pole, 1. / 3); // a^{1/3}
421 double b_a23 = (U - E) / pole;
422 double arg = a13 * (x + b_a23);
423 // std::cerr<<"\nx = "<<x<<" x_pocz = "<<x_pocz<<" pole = "<<pole<<" a13 = "<<a13<<" b_a23 = "<<b_a23<<" arg_Ai' =
424 // "<<arg<<" inaczej (w x_pocz)"<<(2*masa_p(E)*(y_pocz - E)/(a13*a13))<<" inaczej (w x_kon)"<<(2*masa_p(E)*(y_kon -
425 // E)/(a13*a13));
426 return a13 * gsl_sf_airy_Ai_deriv(arg, GSL_PREC_DOUBLE); // Na razie nie ma
427}
428/*****************************************************************************/
429 double warstwa::Ai_prim_skala(double x, double E) const
430{
431 if(y_kon == y_pocz)
432 {
433 Error err; // MD
434 err << "Zla funkcja!\n";
435 throw err;
436 }
437 double U = y_pocz - pole * x_pocz;
438 double a13 = (pole > 0) ? pow(2 * masa_p(E) * pole, 1. / 3) : -pow(-2 * masa_p(E) * pole, 1. / 3); // a^{1/3}
439 double b_a23 = (U - E) / pole;
440 double arg = a13 * (x + b_a23);
441 // std::cerr<<"\nx = "<<x<<" x_pocz = "<<x_pocz<<" pole = "<<pole<<" a13 = "<<a13<<" b_a23 = "<<b_a23<<" arg_Ai' =
442 // "<<arg<<" inaczej (w x_pocz)"<<(2*masa_p(E)*(y_pocz - E)/(a13*a13))<<" inaczej (w x_kon)"<<(2*masa_p(E)*(y_kon -
443 // E)/(a13*a13));
444 double arg_pocz = a13 * (x_pocz + b_a23);
445 double przeskal = 1.0;
446 if(arg > 0 && arg_pocz > 0)
447 {
448 przeskal = exp(-2 * (pow(arg, 3. / 2) - pow(arg_pocz, 3. / 2)) / 3);
449 }
450 else
451 {
452 if(arg > 0)
453 {
454 przeskal = exp(-2 * pow(arg, 3. / 2) / 3);
455 }
456 if(arg_pocz > 0)
457 {
458 przeskal = exp(2 * pow(arg_pocz, 3. / 2) / 3);
459 }
460 }
461 return a13 * gsl_sf_airy_Ai_deriv_scaled(arg, GSL_PREC_DOUBLE) * przeskal; // Na razie nie ma
462}
463/*****************************************************************************/
464 double warstwa::Bi(double x, double E) const
465{
466 if(y_kon == y_pocz)
467 {
468 Error err; // MD
469 err << "Zla funkcja!\n";
470 throw err;
471 }
472 // rownanie: -f''(x) + (b + ax)f(x) = 0
473 // a = 2m*pole/h^2
474 // b = 2m(U - E)/h^2
475 // rozw f(x) = Ai( (ax+b)/a^{2/3} ) = Ai( a^{1/3} (x + b/a^{2/3}) )
476 double U = y_pocz - pole * x_pocz;
477 double a13 = (pole > 0) ? pow(2 * masa_p(E) * pole, 1. / 3) : -pow(-2 * masa_p(E) * pole, 1. / 3); // a^{1/3}
478 double b_a23 = (U - E) / pole;
479 double arg = a13 * (x + b_a23);
480 // std::cerr<<"\narg_Bi = "<<arg;
481 return gsl_sf_airy_Bi(arg, GSL_PREC_DOUBLE); // Na razie nie ma
482}
483/*****************************************************************************/
484 double warstwa::Bi_skala(double x, double E) const
485{
486 if(y_kon == y_pocz)
487 {
488 Error err; // MD
489 err << "Zla funkcja!\n";
490 throw err;
491 }
492 // rownanie: -f''(x) + (b + ax)f(x) = 0
493 // a = 2m*pole/h^2
494 // b = 2m(U - E)/h^2
495 // rozw f(x) = Ai( (ax+b)/a^{2/3} ) = Ai( a^{1/3} (x + b/a^{2/3}) )
496 double U = y_pocz - pole * x_pocz;
497 double a13 = (pole > 0) ? pow(2 * masa_p(E) * pole, 1. / 3) : -pow(-2 * masa_p(E) * pole, 1. / 3); // a^{1/3}
498 double b_a23 = (U - E) / pole;
499 double arg = a13 * (x + b_a23);
500 // std::cerr<<"\narg_Bi = "<<arg;
501 double arg_pocz = a13 * (x_pocz + b_a23);
502 double przeskal = 1.0;
503 if(arg > 0 && arg_pocz > 0)
504 {
505 przeskal = exp(-2 * (pow(arg, 3. / 2) - pow(arg_pocz, 3. / 2)) / 3);
506 }
507 else
508 {
509 if(arg > 0)
510 {
511 przeskal = exp(-2 * pow(arg, 3. / 2) / 3);
512 }
513 if(arg_pocz > 0)
514 {
515 przeskal = exp(2 * pow(arg_pocz, 3. / 2) / 3);
516 }
517 }
518 return gsl_sf_airy_Bi_scaled(arg, GSL_PREC_DOUBLE) / przeskal; // Na razie nie ma
519}
520/*****************************************************************************/
521 double warstwa::Bi_prim(double x, double E) const
522{
523 if(y_kon == y_pocz)
524 {
525 Error err; // MD
526 err << "Zla funkcja!\n";
527 throw err;
528 }
529 double U = y_pocz - pole * x_pocz;
530 double a13 = (pole > 0) ? pow(2 * masa_p(E) * pole, 1. / 3) : -pow(-2 * masa_p(E) * pole, 1. / 3); // a^{1/3}
531 double b_a23 = (U - E) / pole;
532 double arg = a13 * (x + b_a23);
533 // std::cerr<<"\narg_Bi' = "<<arg;
534 return a13 * gsl_sf_airy_Bi_deriv(arg, GSL_PREC_DOUBLE); // Na razie nie ma
535}
536/*****************************************************************************/
537 double warstwa::Bi_prim_skala(double x, double E) const
538{
539 if(y_kon == y_pocz)
540 {
541 Error err; // MD
542 err << "Zla funkcja!\n";
543 throw err;
544 }
545 double U = y_pocz - pole * x_pocz;
546 double a13 = (pole > 0) ? pow(2 * masa_p(E) * pole, 1. / 3) : -pow(-2 * masa_p(E) * pole, 1. / 3); // a^{1/3}
547 double b_a23 = (U - E) / pole;
548 double arg = a13 * (x + b_a23);
549 // std::cerr<<"\narg_Bi' = "<<arg;
550 double arg_pocz = a13 * (x_pocz + b_a23);
551 double przeskal = 1.0;
552 if(arg > 0 && arg_pocz > 0)
553 {
554 przeskal = exp(-2 * (pow(arg, 3. / 2) - pow(arg_pocz, 3. / 2)) / 3);
555 }
556 else
557 {
558 if(arg > 0)
559 {
560 przeskal = exp(-2 * pow(arg, 3. / 2) / 3);
561 }
562 if(arg_pocz > 0)
563 {
564 przeskal = exp(2 * pow(arg_pocz, 3. / 2) / 3);
565 }
566 }
567 return a13 * gsl_sf_airy_Bi_deriv_scaled(arg, GSL_PREC_DOUBLE) / przeskal; // Na razie nie ma
568}
569/*****************************************************************************/
570 double warstwa::airy_kwadr_pierwotna(double x, double E, double A, double B) const
571{
572 if(y_kon == y_pocz)
573 {
574 Error err; // MD
575 err << "Zla funkcja!\n";
576 throw err;
577 }
578 double U = y_pocz - pole * x_pocz;
579 double b_a23 = (U - E) / pole;
580 double a = 2 * masa_p(E) * pole;
581 double f = funkcjafal(x, E, A, B);
582 double fp = funkcjafal_prim(x, E, A, B);
583 return (x + b_a23) * f * f - fp * fp / a;
584}
585/*****************************************************************************/
586 double warstwa::k_kwadr(double E) const // Zwraca k^2, ujemne dla energii spod bariery (- kp^2)
587{
588 double wartosc;
589 if(pole != 0)
590 {
591 Error err; // MD
592 err << "Jesze nie ma airych!\n";
593 throw err;
594 }
595 else
596 {
597 wartosc = 2 * masa_p(E) * (E - y_pocz);
598 }
599 return wartosc;
600}
601/*****************************************************************************/
602 int warstwa::zera_ffal(double E, double A, double B, double sasiad_z_lewej, double sasiad_z_prawej)
603 const // wartosci sasiadow po to, zeby uniknac klopotow, kiedy zero wypada na laczeniu
604{
605 int tylezer = 0;
606 double wart_kon =
607 (funkcjafal(x_kon, E, A, B) + sasiad_z_prawej) / 2; // Usrednienie dla unikniecia klopotow z zerami na laczeniu,
608 // gdzie malutka nieciaglosc moze generowac zmiany znakow
609 double wart_pocz = (funkcjafal(x_pocz, E, A, B) + sasiad_z_lewej) / 2;
610 double iloczyn = wart_pocz * wart_kon;
611 // std::cerr<<"\nwart na koncach: "<<funkcjafal(x_pocz, E, A, B)<<", "<<funkcjafal(x_kon, E, A, B);
612 // std::cerr<<"\npo usrednieniu: "<<wart_pocz<<", "<<wart_kon;
613 if(pole != 0)
614 {
615 double a13 = (pole > 0) ? pow(2 * masa_p(E) * pole, 1. / 3) : -pow(-2 * masa_p(E) * pole, 1. / 3); // a^{1/3}
616 double U = y_pocz - pole * x_pocz;
617 double b_a23 = (U - E) / pole;
618 double arg1, arg2, argl, argp, x1, x2, xlew, xpra;
619 int nrza, nrprzed; // nrza do argp, nrprzed do argl
620 arg1 = a13 * (x_pocz + b_a23);
621 arg2 = a13 * (x_kon + b_a23);
622 // argl = std::min(arg1, arg2); // LUKASZ dodalem komentarz
623 // argp = std::max(arg1, arg2); // LUKASZ dodalem komentarz
624 argl = min(arg1, arg2); // LUKASZ moje funkcje
625 argp = max(arg1, arg2); // LUKASZ moje funkcje
626 nrza = 1;
627 double z1 = -1.174; // oszacowanie pierwszego zera B1
628 double dz = -2.098; // oszacowanie odstepu miedzy perwszymi dwoma zerami
629 // nrprzed=1;
630 // nrprzed = floor((argl-z1)/dz + 1); // oszacowanie z dolu numeru miejsca zerowego
631 nrprzed = floor((argp - z1) / dz + 1);
632 nrprzed = (nrprzed >= 1) ? nrprzed : 1;
633 int tymcz = 0;
634 double ntezero = gsl_sf_airy_zero_Bi(nrprzed);
635 // std::cerr<<"\nU = "<<U<<" a13 = "<<a13<<" b_a23 = "<<b_a23<<" argl = "<<argl<<" argp = "<<argp<<" ntezero
636 // = "<<ntezero<<" nrprzed = "<<nrprzed;
637 double brak; // oszacowanie z dolu braku
638 long licznik = 0;
639 // while(ntezero>=argl)
640 while(ntezero >= argp)
641 {
642 if(nrprzed > 2)
643 {
644 dz = ntezero - gsl_sf_airy_zero_Bi(nrprzed - 1);
645 brak = (argp - ntezero) / dz;
646 if(brak > 2.) // jesli jeszcze daleko
647 {
648 nrprzed = nrprzed + floor(brak);
649 }
650 else
651 nrprzed++;
652 }
653 else
654 nrprzed++;
655 ntezero = gsl_sf_airy_zero_Bi(nrprzed);
656 licznik++;
657 // std::cerr<<"\nnrprzed = "<<nrprzed<<" ntezero = "<<ntezero;
658 }
659 // std::cerr<<"\tnrprzed kon "<<nrprzed<<" po "<<licznik<<" dodawaniach\n";
660 nrza = nrprzed;
661 nrprzed--;
662 // while(gsl_sf_airy_zero_Bi(nrza)>=argp)
663 while(gsl_sf_airy_zero_Bi(nrza) >= argl)
664 {
665 nrza++;
666 // std::cerr<<"\nnrza = "<<nrza<<" ntezero = "<<gsl_sf_airy_zero_Bi(nrza);
667 }
668 // std::cerr<<"\nnrprzed = "<<nrprzed<<" nrza = "<<nrza;
669 /*
670 std::cerr<<"\nnrprzed = "<<nrprzed<<" nrza = "<<nrza;
671 x1 = b_a23 - gsl_sf_airy_zero_Bi(nrprzed+1)/a13; // polozenia skrajnych zer Ai w studni
672 x2 = b_a23 - gsl_sf_airy_zero_Bi(nrza-1)/a13;
673 xlew = std::min(x1, x2);
674 xpra = std::max(x1, x2);
675 std::cerr<<"\txlew="<<struktura::dlugosc_na_A(xlew)<<" xpra="<<struktura::dlugosc_na_A(xpra);
676 tylko do testow tutaj */
677
678 if(nrza - nrprzed >= 2)
679 {
680 tymcz = nrza - nrprzed - 2;
681 x1 = -b_a23 + gsl_sf_airy_zero_Bi(nrprzed + 1) / a13; // polozenia skrajnych zer Ai w studni
682 x2 = -b_a23 + gsl_sf_airy_zero_Bi(nrza - 1) / a13;
683 // xlew = std::min(x1, x2); // LUKASZ dodalem komentarz
684 // xpra = std::max(x1, x2); // LUKASZ dodalem komentarz
685 xlew = min(x1, x2); // LUKASZ moje funkcje
686 xpra = max(x1, x2); // LUKASZ moje funkcje
687 // std::cerr<<"\n xlew="<<struktura::dlugosc_na_A(xlew)<<" xpra="<<struktura::dlugosc_na_A(xpra);
688 // std::cerr<<"\n A "<<funkcja_z_polem_do_oo(struk.punkty[i],E,funkcja,struk)<<"
689 // "<<funkcja_z_polem_do_oo(xl,E,funkcja,struk);
690 if(wart_pocz * funkcjafal(xlew, E, A, B) < 0)
691 // if(funkcja_z_polem_do_oo(struk.punkty[i],E,funkcja,struk)*funkcja_z_polem_do_oo(xl,E,funkcja,struk)<0)
692 tymcz++;
693 if(wart_kon * funkcjafal(xpra, E, A, B) < 0)
694 tymcz++;
695 }
696 else
697 {
698 // std::cerr<<"\n C "<<funkcja_z_polem_do_oo(struk.punkty[i+1],E,funkcja,struk)<<"
699 // "<<funkcja_z_polem_do_oo(struk.punkty[i],E,funkcja,struk);
700 if(iloczyn < 0)
701 tymcz = 1;
702 }
703 tylezer = tymcz;
704 // Error err; // MD
705 // err<<"Jeszcze nie ma zer Airy'ego!\n";
706 // throw err;
707 }
708 else
709 {
710 if(E >= y_pocz)
711 {
712 double k = sqrt(2 * masa_p(E) * (E - y_pocz));
713 tylezer = int(k * (x_kon - x_pocz) / M_PI);
714 if(tylezer % 2 == 0)
715 {
716 if(iloczyn < 0)
717 {
718 tylezer++;
719 }
720 }
721 else
722 {
723 if(iloczyn > 0)
724 {
725 tylezer++;
726 }
727 }
728 }
729 else
730 {
731 if(iloczyn < 0)
732 {
733 tylezer++;
734 }
735 }
736 }
737 // std::cerr<<"\nE = "<<E<<"\tiloczyn = "<<iloczyn<<"\t zer jest "<<tylezer;
738 // std::cerr<<"\nE = "<<E<<"\t zer jest "<<tylezer;
739 return tylezer;
740}
741/*****************************************************************************/
742 int warstwa::zera_ffal(double E, double A, double B) const
743{
744 int tylezer = 0;
745 double wart_kon = funkcjafal(x_kon, E, A, B);
746 double iloczyn = funkcjafal(x_pocz, E, A, B) * wart_kon;
747 // std::cerr<<"\n wart na koncach: "<<funkcjafal(x_pocz, E, A, B)<<", "<<funkcjafal(x_kon, E, A, B);
748 if(pole != 0)
749 {
750 double a13 = (pole > 0) ? pow(2 * masa_p(E) * pole, 1. / 3) : -pow(-2 * masa_p(E) * pole, 1. / 3); // a^{1/3}
751 double U = y_pocz - pole * x_pocz;
752 double b_a23 = (U - E) / pole;
753 double arg1, arg2, argl, argp, x1, x2, xlew, xpra;
754 int nrza, nrprzed; // nrza do argp, nrprzed do argl
755 arg1 = a13 * (x_pocz + b_a23);
756 arg2 = a13 * (x_kon + b_a23);
757 // argl = std::min(arg1, arg2); // LUKASZ dodalem komentarz
758 // argp = std::max(arg1, arg2); // LUKASZ dodalem komentarz
759 argl = min(arg1, arg2); // LUKASZ moje funkcje
760 argp = max(arg1, arg2); // LUKASZ moje funkcje
761 nrza = 1;
762 double z1 = -1.174; // oszacowanie pierwszego zera B1
763 double dz = -2.098; // oszacowanie odstepu miedzy perwszymi dwoma zerami
764 // nrprzed=1;
765 // nrprzed = floor((argl-z1)/dz + 1); // oszacowanie z dolu numeru miejsca zerowego
766 nrprzed = floor((argp - z1) / dz + 1);
767 nrprzed = (nrprzed >= 1) ? nrprzed : 1;
768 int tymcz = 0;
769 double ntezero = gsl_sf_airy_zero_Bi(nrprzed);
770 // std::cerr<<"\nU = "<<U<<" a13 = "<<a13<<" b_a23 = "<<b_a23<<" argl = "<<argl<<" argp = "<<argp<<" ntezero
771 // = "<<ntezero<<" nrprzed = "<<nrprzed;
772 double brak; // oszacowanie z dolu braku
773 long licznik = 0;
774 // while(ntezero>=argl)
775 while(ntezero >= argp)
776 {
777 if(nrprzed > 2)
778 {
779 dz = ntezero - gsl_sf_airy_zero_Bi(nrprzed - 1);
780 brak = (argp - ntezero) / dz;
781 if(brak > 2.) // jesli jeszcze daleko
782 {
783 nrprzed = nrprzed + floor(brak);
784 }
785 else
786 nrprzed++;
787 }
788 else
789 nrprzed++;
790 ntezero = gsl_sf_airy_zero_Bi(nrprzed);
791 licznik++;
792 // std::cerr<<"\nnrprzed = "<<nrprzed<<" ntezero = "<<ntezero;
793 }
794 // std::cerr<<"\tnrprzed kon "<<nrprzed<<" po "<<licznik<<" dodawaniach\n";
795 nrza = nrprzed;
796 nrprzed--;
797 // while(gsl_sf_airy_zero_Bi(nrza)>=argp)
798 while(gsl_sf_airy_zero_Bi(nrza) >= argl)
799 {
800 nrza++;
801 // std::cerr<<"\nnrza = "<<nrza<<" ntezero = "<<gsl_sf_airy_zero_Bi(nrza);
802 }
803 // std::cerr<<"\nnrprzed = "<<nrprzed<<" nrza = "<<nrza;
804 /*
805 std::cerr<<"\nnrprzed = "<<nrprzed<<" nrza = "<<nrza;
806 x1 = b_a23 - gsl_sf_airy_zero_Bi(nrprzed+1)/a13; // polozenia skrajnych zer Ai w studni
807 x2 = b_a23 - gsl_sf_airy_zero_Bi(nrza-1)/a13;
808 xlew = std::min(x1, x2);
809 xpra = std::max(x1, x2);
810 std::cerr<<"\txlew="<<struktura::dlugosc_na_A(xlew)<<" xpra="<<struktura::dlugosc_na_A(xpra);
811 tylko do testow tutaj */
812
813 if(nrza - nrprzed >= 2)
814 {
815 tymcz = nrza - nrprzed - 2;
816 x1 = -b_a23 + gsl_sf_airy_zero_Bi(nrprzed + 1) / a13; // polozenia skrajnych zer Ai w studni
817 x2 = -b_a23 + gsl_sf_airy_zero_Bi(nrza - 1) / a13;
818 // xlew = std::min(x1, x2); // LUKASZ dodalem komentarz
819 // xpra = std::max(x1, x2); // LUKASZ dodalem komentarz
820 xlew = min(x1, x2); // LUKASZ moje funkcje
821 xpra = max(x1, x2); // LUKASZ moje funkcje
822#ifdef LOGUJ // MD
823 std::cerr << "\n xlew=" << struktura::dlugosc_na_A(xlew) << " xpra=" << struktura::dlugosc_na_A(xpra);
824 // std::cerr << "\n A " << funkcja_z_polem_do_oo(struk.punkty[i], E, funkcja, struk) << " " <<
825 // funkcja_z_polem_do_oo(xl, E, funkcja, struk);
826#endif
827 if(funkcjafal(x_pocz, E, A, B) * funkcjafal(xlew, E, A, B) < 0)
828 // if(funkcja_z_polem_do_oo(struk.punkty[i],E,funkcja,struk)*funkcja_z_polem_do_oo(xl,E,funkcja,struk)<0)
829 tymcz++;
830 if(wart_kon * funkcjafal(xpra, E, A, B) < 0)
831 tymcz++;
832 }
833 else
834 {
835 // std::cerr<<"\n C "<<funkcja_z_polem_do_oo(struk.punkty[i+1],E,funkcja,struk)<<"
836 // "<<funkcja_z_polem_do_oo(struk.punkty[i],E,funkcja,struk);
837 if(funkcjafal(x_pocz, E, A, B) * wart_kon < 0)
838 tymcz = 1;
839 }
840 tylezer = tymcz;
841 // Error err; // MD
842 // err<<"Jeszcze nie ma zer Airy'ego!\n";
843 // throw err;
844 }
845 else
846 {
847 if(E >= y_pocz)
848 {
849 double k = sqrt(2 * masa_p(E) * (E - y_pocz));
850 tylezer = int(k * (x_kon - x_pocz) / M_PI);
851 if(tylezer % 2 == 0)
852 {
853 if(iloczyn < 0)
854 {
855 tylezer++;
856 }
857 }
858 else
859 {
860 if(iloczyn > 0)
861 {
862 tylezer++;
863 }
864 }
865 }
866 else
867 {
868 if(iloczyn < 0)
869 {
870 tylezer++;
871 }
872 }
873 }
874#ifdef LOGUJ // MD
875 std::cerr << "\nE = " << E << "\tiloczyn = " << iloczyn << "\t zer jest " << tylezer;
876#endif
877 return tylezer;
878}
879/*****************************************************************************/
880 double warstwa::funkcjafal(double x, double E, double A, double B) const { return A * ffala(x, E) + B * ffalb(x, E); }
881/*****************************************************************************/
882 double warstwa::funkcjafal_prim(double x, double E, double A, double B) const
883{
884 return A * ffala_prim(x, E) + B * ffalb_prim(x, E);
885}
886/*****************************************************************************/
887// Nowe dopiski do wersji 1.5
888/*****************************************************************************/
889 stan::stan(double E, std::vector<parad> & W, int lz)
890{
891 poziom = E;
892 int M = W.size();
893 wspolczynniki.resize(2*(M-2) + 2);
894 wspolczynniki[0] = W[0].second;
895 for(int i = 1; i <= M-2; i++)
896 {
897 wspolczynniki[2*i-1] = W[i].first;
898 wspolczynniki[2*i] = W[i].second;
899 }
900 wspolczynniki[2*M-3] = W[M-1].first;
901 liczba_zer = lz;
902 prawdopodobienstwa.reserve(M/2 +1);
903}
904
905
906 parad warstwa::AB_z_wartp(double w, double wp, double E) const // liczy ampitudy A i B na podstawie wartości w początku warstwy. wp to pochodna przez masę
907{
908 double masa = masa_p(E);
909 double mian = (ffala_prim(x_pocz, E)*ffalb(x_pocz, E) - ffalb_prim(x_pocz, E)*ffala(x_pocz, E));
910 double A = (-w*ffalb_prim(x_pocz, E) + masa*wp*ffalb(x_pocz, E)) / mian;
911 double B = - (-w*ffala_prim(x_pocz, E) + masa*wp*ffala(x_pocz, E)) / mian;
912 // std::cerr<<"m_p i masa_p: "<<m_p<<", "<<masa_p(E)<<"\n";
913 return {A, B};
914}
915
916 parad warstwa_skraj::AB_z_wartp(double w, double wp, double E) const // liczy ampitudy A i B na podstawie wartości w początku warstwy. wp to pochodna przez masę
917{
918 double masa = masa_p;
919 double mian = (expa_prim(iks, E)*expb(iks, E) - expb_prim(iks, E)*expa(iks, E));
920 double A = (-w*expb_prim(iks, E) + masa*wp*expb(iks, E)) / mian;
921 double B = - (-w*expa_prim(iks, E) + masa*wp*expa(iks, E)) / mian;
922 // std::cerr<<"skrajna:\niks = "<<iks<<", x_pocz = "<<x_pocz<<"\n";
923 // std::cerr<<"skrajna:\nexpa = "<<expa(iks, E)<<", expb = "<<expb(iks, E)<<"\n";
924 // std::cerr<<"skrajna:\nexpa' = "<<expa_prim(iks, E)<<", expb' = "<<expb_prim(iks, E)<<"\n";
925 return {A, B};
926}
927
928 parad struktura::sklejanie_od_lewej(double E)
929{
930 double A, B, w, wp;
931 warstwa war = kawalki[0]; // bo nie ma konstruktora domyślnego
932 // std::cerr<<"lewa\n";
933 w = lewa.funkcjafal(lewa.iks, E, 1.0); // wartosci dla wpółczynnika 1
934 wp = lewa.funkcjafal_prim(lewa.iks, E, 1.0)/lewa.masa_p;
935 // std::cerr<<"w, wp: "<<w<<", "<<wp<<"\n";
936 parad AB;
937 for(int i = 0; i <= (int)kawalki.size() - 1; i++)
938 {
939 // std::cerr<<i<<". warstwa\n";
940 war = kawalki[i];
941 AB = war.AB_z_wartp(w, wp, E);
942 A = AB.first;
943 B = AB.second;
944 // std::cerr<<"AB: "<<A<<", "<<B<<"\n";
945 w = war.funkcjafal(war.x_kon, E, A, B);
946 wp = war.funkcjafal_prim(war.x_kon, E, A, B)/war.masa_p(E);
947 // std::cerr<<"lewe w, wp: "<<war.funkcjafal(war.x_pocz, E, A, B)<<", "<<war.funkcjafal_prim(war.x_pocz, E, A, B)/war.masa_p(E)<<"\n";
948 // std::cerr<<"w, wp: "<<w<<", "<<wp<<"\n";
949 }
950 // std::cerr<<"prawa\n";
951 AB = prawa.AB_z_wartp(w, wp, E);
952 A = AB.first;
953 B = AB.second;
954
955 // std::cerr<<"lewe w, wp: "<<prawa.warstwa::funkcjafal(prawa.x_pocz, E, A, B)<<", "<<prawa.warstwa::funkcjafal_prim(prawa.x_pocz, E, A, B)/prawa.masa_p<<"\n";
956 return AB;
957}
958
959std::vector<parad> struktura::wsp_sklejanie_od_lewej(double E)
960{
961 double A, B, w, wp;
962 std::vector<parad> wspolczynniki;
963 warstwa war = kawalki[0]; // bo nie ma konstruktora domyślnego
964 w = lewa.funkcjafal(lewa.iks, E, 1.0); // wartosci dla wpółczynnika 1
965 wp = lewa.funkcjafal_prim(lewa.iks, E, 1.0)/lewa.masa_p;
966 parad AB;
967 wspolczynniki.push_back(parad(0.0, 1.0));
968 for(int i = 0; i <= (int)kawalki.size() - 1; i++)
969 {
970 war = kawalki[i];
971 AB = war.AB_z_wartp(w, wp, E);
972 wspolczynniki.push_back(AB);
973 A = AB.first;
974 B = AB.second;
975 w = war.funkcjafal(war.x_kon, E, A, B);
976 wp = war.funkcjafal_prim(war.x_kon, E, A, B)/war.masa_p(E);
977 }
978 AB = prawa.AB_z_wartp(w, wp, E);
979 wspolczynniki.push_back(AB);
980 return wspolczynniki;
981}
982
983 int struktura::ilezer_ffal(double E, std::vector<parad> & W)
984{
985 int N = kawalki.size() + 2; //liczba warstw
986 int sumazer = 0;
987 double A, B;//, Al, Bl, Ap, Bp;
988 int pierwsza = -1;
989 do
990 {
991 pierwsza++;
992 }while(pierwsza <= N-3 && (kawalki[pierwsza].y_pocz > E && kawalki[pierwsza].y_kon > E) );
993 int ostatnia = N-2;
994 do
995 {
996 ostatnia--;
997 }while(ostatnia >= 0 && (kawalki[ostatnia].y_pocz > E && kawalki[ostatnia].y_kon > E) );
998 // double sasiadl, sasiadp; // wartość ffal na lewym brzegu sąsiada z prawej (ta, która powinna być wspólna do obu)
999 if(ostatnia == pierwsza) // tylko jedna podejrzana warstwa, nie trzeba się sąsiadami przejmować
1000 {
1001 A = W[pierwsza+1].first;
1002 B = W[pierwsza+1].second;
1003 sumazer += kawalki[pierwsza].zera_ffal(E, A, B);
1004 }
1005 else
1006 {
1007 int j = pierwsza;
1008 A = W[j+1].first;
1009 B = W[j+1].second;
1010
1011 /*
1012 A = V[2*j+1][V.dim2() - 1];
1013 B = V[2*j+2][V.dim2() - 1];
1014 Ap = V[2*(j+1)+1][V.dim2() - 1];
1015 Bp = V[2*(j+1)+2][V.dim2() - 1];
1016 sasiadp = kawalki[j+1].funkcjafal(kawalki[j+1].x_pocz, E, Ap, Bp);
1017 sasiadl = kawalki[j].funkcjafal(kawalki[j].x_pocz, E, A, B); // po lewej nie ma problemu, więc można podstawić wartość z własnej warstwy
1018 sumazer += kawalki[j].zera_ffal(E, A, B, sasiadl, sasiadp);
1019 */
1020 sumazer += kawalki[j].zera_ffal(E, A, B);
1021 // for(int j = pierwsza + 1; j <= ostatnia - 1; j++)
1022 for(int j = pierwsza + 1; j <= ostatnia; j++)
1023 {
1024 A = W[j+1].first;
1025 B = W[j+1].second;
1026 sumazer += kawalki[j].zera_ffal(E, A, B);
1027 /*
1028 Al = V[2*(j-1)+1][V.dim2() - 1];
1029 Bl = V[2*(j-1)+2][V.dim2() - 1];
1030 A = V[2*j+1][V.dim2() - 1];
1031 B = V[2*j+2][V.dim2() - 1];
1032 Ap = V[2*(j+1)+1][V.dim2() - 1];
1033 Bp = V[2*(j+1)+2][V.dim2() - 1];
1034 sasiadl = kawalki[j-1].funkcjafal(kawalki[j-1].x_kon, E, Al, Bl);
1035 sasiadp = kawalki[j+1].funkcjafal(kawalki[j+1].x_pocz, E, Ap, Bp);
1036 sumazer += kawalki[j].zera_ffal(E, A, B, sasiadl, sasiadp);
1037 */
1038 }
1039 /*
1040 j = ostatnia;
1041 A = V[2*j+1][V.dim2() - 1];
1042 B = V[2*j+2][V.dim2() - 1];
1043 Al = V[2*(j-1)+1][V.dim2() - 1];
1044 Bl = V[2*(j-1)+2][V.dim2() - 1];
1045 sasiadp = kawalki[j].funkcjafal(kawalki[j].x_kon, E, A, B); // po prawej nie ma problemu, więc można podstawić wartość z własnej warstwy
1046 sasiadl = kawalki[j-1].funkcjafal(kawalki[j-1].x_kon, E, Al, Bl);
1047 sumazer += kawalki[j].zera_ffal(E, A, B, sasiadl, sasiadp); // W ostatniej warswie nie może być zera na łączeniu, więc nie ma problemu
1048 */
1049 }
1050 return sumazer;
1051}
1052
1053
1054 double struktura::blad_sklejania(double E)
1055{
1056 return sklejanie_od_lewej(E).second; // drugi współczynnk powinien być 0
1057}
1058
1059 double struktura::poprawianie_poziomu(double E, double DE)
1060{
1061 double B0 = blad_sklejania(E);
1062 double E1 = E + DE;
1063 double B1 = blad_sklejania(E1);
1064 // std::cerr<<"poprawianie. B0 i B1: "<<B0<<", "<<B1<<"\n";
1065 if(B0*B1 > 0 && abs(B0) < abs(B1)) // idziemy w złym kierunku
1066 {
1067 DE = - DE; //odwracamy i zakadamy, że teraz będzie dobrze. ryzykownie
1068 E1 = E + DE;
1069 B1 = blad_sklejania(E1);
1070 // std::cerr<<"zmiana znaku. B0 i B1: "<<B0<<", "<<B1<<"\n";
1071 }
1072
1073 while(B0*B1 > 0 && abs(B0) > abs(B1))
1074 {
1075 E1 += DE;
1076 B1 = blad_sklejania(E1);
1077 // std::cerr<<"drążymy. B0 i B1: "<<B0<<", "<<B1<<"\n";
1078 }
1079 // std::cerr<<"poprawianie. B0 i B1: "<<B0<<", "<<B1<<"\n";
1080 double (struktura::*fun)(double) = & struktura::blad_sklejania;
1081 double wyn = bisekcja(fun, (std::min)(E, E1), (std::max)(E, E1), 1e-15); // trzeba dokładniej niż domyślnie
1082 return wyn;
1083}
1084
1085 void struktura::szukanie_poziomow_zpoprawka(double Ek, double rozdz)
1086{
1087 std::list<stan> listast; // lista znalezionych stanów
1088 std::list<stan>::iterator itostd, itnastd; // najwyższy stan nad któ©ym wszystko jest znalezione, następny stan, nad którym trzeba jeszcze szukać
1089 std::list<stan>::iterator iter; //pomocniczy
1090 double E0 = dol;
1091 if(Ek <= E0)
1092 {
1093 std::cerr<<"Zły zakres energii!\n";
1094 abort();
1095 }
1096 int M = 2*(kawalki.size() + 2) - 2;
1097 double wartakt;
1098 double (struktura::*fun)(double) = & struktura::czyosobliwa;
1099 std::vector<double> trojka;
1100 std::clog<<"W szukaniu, E_pocz = "<<Ek<<"\n";
1101 double wartpop = czyosobliwa(Ek);
1102 std::clog<<"Pierwsza wartosc = "<<wartpop<<"\n";
1103 A2D V(M, M);
1104 //1.5
1105 std::vector<parad> wspolcz;
1106 //
1107 int ilepoziomow;
1108 int liczbazer;
1109 double E, Epop;
1110 // double B0, B1;
1111 double mnozdorozdz = 1e3; //monnik zmiennej rozdz do kroku w szukaniu górnego stanu
1112 if(wartpop == 0)
1113 {
1114 Ek -= rozdz;
1115 wartpop = czyosobliwa(Ek);
1116 }
1117 Epop = Ek;
1118 E = Ek - rozdz*mnozdorozdz;
1119 wartakt = czyosobliwa(E);
1120 while( (E > E0) && (wartpop*wartakt > 0) )
1121 {
1122 wartpop = wartakt;
1123 Epop = E;
1124 E -= rozdz*mnozdorozdz; // żeby nie czekać wieki
1125 wartakt = czyosobliwa(E);
1126 }
1127 if(E <= E0)
1128 {
1129 std::cerr<<"Nie znalazł nawet jednego poziomu!\n";
1130 abort();
1131 }
1132 std::clog<<"Bisekcja z krancami "<<E<<" i "<<Epop<<"\n";
1133 double Eost = bisekcja(fun, E, Epop); // wstępna energia najwyższego poziomu
1134
1135 // dopiski 1.5:
1136 double nowaE = poprawianie_poziomu(Eost, 1e-11);
1137 wspolcz = wsp_sklejanie_od_lewej(nowaE);
1138 Eost = nowaE;
1139
1140 liczbazer = ilezer_ffal(Eost, wspolcz);
1141 ilepoziomow = liczbazer + 1; // Tyle ma znaleźć poziomów
1142 stan nowy = stan(Eost, wspolcz, ilepoziomow - 1);
1143
1144 listast.push_back(nowy);
1145 itostd = listast.end();
1146 --itostd; // iterator na ostatni element
1147 stan stymcz;
1148 std::clog<<" Eost = "<<Eost<<" ilepoziomow = "<<ilepoziomow<<"\n";
1149 punkt ost, ostdob, pierw;
1150 itnastd = listast.begin();
1151 // while(listast.size() < ilepoziomow)
1152 // while(itostd != listast.begin())
1153 while((int)listast.size() < ilepoziomow)
1154 {
1155 // licz = 0;
1156 ostdob = *itostd;
1157 E = ostdob.en - rozdz;
1158 ost = punkt(E, czyosobliwa(E));
1159 // E = (nast_dobre >= 0)?(rozwiazania[nast_dobre].poziom + rozdz):(E0 + rozdz);
1160 E = (itnastd != itostd)?(itnastd->poziom + rozdz):(E0 + rozdz);
1161 pierw = punkt(E, czyosobliwa(E));
1162
1163 if(ost.wart*pierw.wart > 0)
1164 {
1165 std::clog<<"Zagęszczanie z pierw = "<<pierw.en<<" i ost = "<<ost.en<<"\n";
1166 trojka = zageszczanie(pierw, ost);
1167 if(!trojka.empty())
1168 {
1169 iter = itostd;
1170 for(int i = 1; i >= 0; i--)// bo są dwa zera
1171 {
1172 E = bisekcja(fun, trojka[i], trojka[i+1]);
1173 nowaE = poprawianie_poziomu(E, 1e-11);
1174 E = nowaE;
1175 wspolcz = wsp_sklejanie_od_lewej(E);
1176 liczbazer = ilezer_ffal(E, wspolcz);
1177 nowy = stan(E, wspolcz, liczbazer);
1178 std::clog<<"E = "<<E<<" zer jest "<<liczbazer<<"\n";
1179 listast.insert(iter, nowy);
1180 --iter;
1181 }
1182 }
1183 else
1184 {
1185 itostd = itnastd;
1186 std::clog<<"Nie znalazł a powinien\n";
1187 }
1188
1189 }
1190 else
1191 {
1192 E = bisekcja(fun, pierw.en, ost.en);
1193 nowaE = poprawianie_poziomu(E, 1e-11);
1194 E = nowaE;
1195 wspolcz = wsp_sklejanie_od_lewej(E);
1196 liczbazer = ilezer_ffal(E, wspolcz);
1197 nowy = stan(E, wspolcz, liczbazer);
1198 listast.insert(itostd, nowy);
1199 }
1200 while((itostd != listast.begin()) && (itostd->liczba_zer <= std::prev(itostd)->liczba_zer + 1)) //stany na liście schodzą po jeden w liczbie zer
1201 {
1202 if(itostd->liczba_zer < std::prev(itostd)->liczba_zer + 1)
1203 {
1204 std::clog<<"\nKłopot z miejscami zerowymi: E1 = "<<(itostd->poziom)<<" zer "<<(itostd->liczba_zer)<<"\nE2 = "<<(std::prev(itostd)->poziom)<<" zer "<<(std::prev(itostd)->liczba_zer)<<"\n";
1205 }
1206 --itostd;
1207 }
1208
1209 itnastd = (itostd == listast.begin())?listast.begin():std::prev(itostd);
1210 std::clog<<"ostatnie_dobre = "<<(itostd->liczba_zer)<<" nastepne_dobre = "<<(itnastd->liczba_zer)<<"\n";
1211 }
1212
1213 std::clog<<"\n";
1214 rozwiazania.reserve(listast.size());
1215 iter = listast.begin();
1216 while(iter != listast.end())
1217 {
1218 rozwiazania.push_back(*iter);
1219 std::clog<<"zera = "<<(iter->liczba_zer)<<"\tE = "<<(iter->poziom)<<"\n";
1220 ++iter;
1221 }
1222}
1223
1224 void struktura::szukanie_poziomow_zpoprawka2(double Ek, double rozdz)
1225{
1226 std::list<stan> listast; // lista znalezionych stanów
1227 std::list<stan>::iterator itostd, itnastd; // najwyższy stan nad któ©ym wszystko jest znalezione, następny stan, nad którym trzeba jeszcze szukać
1228 std::list<stan>::iterator iter; //pomocniczy
1229 double E0 = dol;
1230 if(Ek <= E0)
1231 {
1232 std::cerr<<"Zły zakres energii!\n";
1233 abort();
1234 }
1235 int M = 2*(kawalki.size() + 2) - 2;
1236 double wartakt;
1237 double (struktura::*fun)(double) = & struktura::czyosobliwa;
1238 std::vector<double> trojka;
1239 std::clog<<"W szukaniu, E_pocz = "<<Ek<<"\n";
1240 double wartpop = czyosobliwa(Ek);
1241 std::clog<<"Pierwsza wartosc = "<<wartpop<<"\n";
1242 A2D V(M, M);
1243 //1.5
1244 std::vector<parad> wspolcz;
1245 //
1246 int ilepoziomow;
1247 int liczbazer;
1248 double E, Epop;
1249 // double B0, B1;
1250 double mnozdorozdz = 1e3; //monnik zmiennej rozdz do kroku w szukaniu górnego stanu
1251 if(wartpop == 0)
1252 {
1253 Ek -= rozdz;
1254 wartpop = czyosobliwa(Ek);
1255 }
1256 Epop = Ek;
1257 E = Ek - rozdz*mnozdorozdz;
1258 wartakt = czyosobliwa(E);
1259 while( (E > E0) && (wartpop*wartakt > 0) )
1260 {
1261 wartpop = wartakt;
1262 Epop = E;
1263 E -= rozdz*mnozdorozdz; // żeby nie czekać wieki
1264 wartakt = czyosobliwa(E);
1265 }
1266 if(E <= E0)
1267 {
1268 std::cerr<<"Nie znalazł nawet jednego poziomu!\n";
1269 abort();
1270 }
1271 std::clog<<"Bisekcja z krancami "<<E<<" i "<<Epop<<"\n";
1272 double Eost = bisekcja(fun, E, Epop); // wstępna energia najwyższego poziomu
1273 double nowaE = poprawianie_poziomu(Eost, 1e-11);
1274 Eost = nowaE;
1275 // V = A2D(M, M);
1276 liczbazer = ilezer_ffal(Eost, V);
1277
1278 // liczbazer = ilezer_ffal(Eost, wspolcz);
1279 ilepoziomow = liczbazer + 1; // Tyle ma znaleźć poziomów
1280 stan nowy = stan(Eost, V, ilepoziomow - 1);
1281
1282 std::cerr<<" Eost = "<<Eost<<" ilepoziomow = "<<ilepoziomow<<"\n";
1283
1284 listast.push_back(nowy);
1285 itostd = listast.end();
1286 --itostd; // iterator na ostatni element
1287 stan stymcz;
1288 std::clog<<" Eost = "<<Eost<<" ilepoziomow = "<<ilepoziomow<<"\n";
1289 punkt ost, ostdob, pierw;
1290 itnastd = listast.begin();
1291 // while(listast.size() < ilepoziomow)
1292 // while(itostd != listast.begin())
1293 while((int)listast.size() < ilepoziomow)
1294 {
1295 // licz = 0;
1296 ostdob = *itostd;
1297 E = ostdob.en - rozdz;
1298 ost = punkt(E, czyosobliwa(E));
1299 // E = (nast_dobre >= 0)?(rozwiazania[nast_dobre].poziom + rozdz):(E0 + rozdz);
1300 E = (itnastd != itostd)?(itnastd->poziom + rozdz):(E0 + rozdz);
1301 pierw = punkt(E, czyosobliwa(E));
1302
1303 if(ost.wart*pierw.wart > 0)
1304 {
1305 std::clog<<"Zagęszczanie z pierw = "<<pierw.en<<" i ost = "<<ost.en<<"\n";
1306 trojka = zageszczanie(pierw, ost);
1307 if(!trojka.empty())
1308 {
1309 iter = itostd;
1310 for(int i = 1; i >= 0; i--)// bo są dwa zera
1311 {
1312 E = bisekcja(fun, trojka[i], trojka[i+1]);
1313 nowaE = poprawianie_poziomu(E, 1e-11);
1314 E = nowaE;
1315 liczbazer = ilezer_ffal(E, V);
1316 /* wspolcz = wsp_sklejanie_od_lewej(E);
1317 liczbazer = ilezer_ffal(E, wspolcz);*/
1318 nowy = stan(E, V, liczbazer);
1319 std::clog<<"E = "<<E<<" zer jest "<<liczbazer<<"\n";
1320 listast.insert(iter, nowy);
1321 --iter;
1322 }
1323 }
1324 else
1325 {
1326 itostd = itnastd;
1327 std::clog<<"Nie znalazł a powinien\n";
1328 }
1329
1330 }
1331 else
1332 {
1333 E = bisekcja(fun, pierw.en, ost.en);
1334 nowaE = poprawianie_poziomu(E, 1e-11);
1335 E = nowaE;
1336 liczbazer = ilezer_ffal(E, V);
1337 /* wspolcz = wsp_sklejanie_od_lewej(E);
1338 liczbazer = ilezer_ffal(E, wspolcz);*/
1339 nowy = stan(E, V, liczbazer);
1340 listast.insert(itostd, nowy);
1341 }
1342 while((itostd != listast.begin()) && (itostd->liczba_zer <= std::prev(itostd)->liczba_zer + 1)) //stany na liście schodzą po jeden w liczbie zer
1343 {
1344 if(itostd->liczba_zer < std::prev(itostd)->liczba_zer + 1)
1345 {
1346 std::clog<<"\nKłopot z miejscami zerowymi: E1 = "<<(itostd->poziom)<<" zer "<<(itostd->liczba_zer)<<"\nE2 = "<<(std::prev(itostd)->poziom)<<" zer "<<(std::prev(itostd)->liczba_zer)<<"\n";
1347 }
1348 --itostd;
1349 }
1350
1351 itnastd = (itostd == listast.begin())?listast.begin():std::prev(itostd);
1352 std::clog<<"ostatnie_dobre = "<<(itostd->liczba_zer)<<" nastepne_dobre = "<<(itnastd->liczba_zer)<<"\n";
1353 }
1354
1355 std::clog<<"\n";
1356 rozwiazania.reserve(listast.size());
1357 iter = listast.begin();
1358 while(iter != listast.end())
1359 {
1360 rozwiazania.push_back(*iter);
1361 std::clog<<"zera = "<<(iter->liczba_zer)<<"\tE = "<<(iter->poziom)<<"\n";
1362 ++iter;
1363 }
1364}
1365
1366 double obszar_aktywny::przekrycia_schodkowe(double E, int nr_c, int nr_v)
1367{
1368 struktura * el = pasmo_przew[nr_c];
1369 struktura * dziu = pasmo_wal[nr_v];
1370 A2D * m_prz = calki_przekrycia[nr_c][nr_v];
1371 double wynik = 0;
1372
1373 double E0;
1374 std::clog<<"E = "<<E<<"\n";
1375 for(int nrpoz_el = 0; nrpoz_el <= int(el->rozwiazania.size()) - 1; nrpoz_el++)
1376 for(int nrpoz_dziu = 0; nrpoz_dziu <= int(dziu->rozwiazania.size()) - 1; nrpoz_dziu++)
1377 {
1378 // std::cerr<<"\nprzekrycie w "<<nrpoz_el<<", "<<nrpoz_dziu;
1379 // std::cerr<<" = "<<(*m_prz)[nrpoz_el][nrpoz_dziu];
1380 E0 = Egcv[nr_v] - Egcc[nr_c] + el->rozwiazania[nrpoz_el].poziom + dziu->rozwiazania[nrpoz_dziu].poziom;
1381 if( E-E0 >= 0 )
1382 {
1383 std::clog<<"E0 = "<<E0<<" dodaję "<<nrpoz_el<<", "<<nrpoz_dziu<<": "<<((*m_prz)[nrpoz_el][nrpoz_dziu])<<"\n";
1384 wynik += ((*m_prz)[nrpoz_el][nrpoz_dziu]);
1385 }
1386 else
1387 break;
1388 }
1389 std::clog<<"wynik = "<<wynik<<"\n";
1390 return wynik;
1391}
1392
1393 void obszar_aktywny::przekrycia_dopliku(ofstream & plik, int nr_c, int nr_v)
1394{
1395 struktura * el = pasmo_przew[nr_c];
1396 struktura * dziu = pasmo_wal[nr_v];
1397 std::vector<parad> Esum;
1398 // double dno = Egcv[nr_v] - Egcc[nr_c] + el->rozwiazania[0].poziom + dziu->rozwiazania[0].poziom // najnuższa energia przejścia
1399 double wart;
1400 double E0;
1401 for(int nrpoz_el = 0; nrpoz_el <= int(el->rozwiazania.size()) - 1; nrpoz_el++)
1402 for(int nrpoz_dziu = 0; nrpoz_dziu <= int(dziu->rozwiazania.size()) - 1; nrpoz_dziu++)
1403 {
1404 E0 = Egcv[nr_v] - Egcc[nr_c] + el->rozwiazania[nrpoz_el].poziom + dziu->rozwiazania[nrpoz_dziu].poziom;
1405 // Esum.push_back({E0, wartpop});
1406 // plik<<E0<<" "<<wartpop<<"\n";
1407 wart = przekrycia_schodkowe(std::nextafter(E0, E0+1), nr_c, nr_v);
1408 Esum.push_back({E0, wart});
1409 // plik<<E0<<" "<<wart<<"\n"; // nextafter bierze następną liczbę double w kieryunku drugiego argumentu
1410 // wartpop = wart;
1411 }
1412 sort(Esum.begin(), Esum.end(), jaksortpar);
1413 double wartpop = 0;
1414 for(int i = 0; i <= (int)Esum.size() - 1; i++)
1415 {
1416 plik<<Esum[i].first<<" "<<wartpop<<"\n";
1417 plik<<Esum[i].first<<" "<<Esum[i].second<<"\n";
1418 wartpop = Esum[i].second;
1419 }
1420}
1421
1422 bool jaksortpar(parad a, parad b) //sortuje po pierwszej wartości, a jak równe , to po drugiej
1423{
1424 bool wyn = 0;
1425 if(a.first < b.first)
1426 wyn = 1;
1427 if(a.first == b.first && a.second < b.second)
1428 wyn = 1;
1429
1430 return wyn;
1431}
1432
1433
1434//koniec dopisków 1.5
1435
1436/*****************************************************************************/
1437 double warstwa::norma_kwadr(double E, double A, double B) const
1438{
1439 double wartosc;
1440 if(pole != 0)
1441 {
1442 wartosc = 1.; // chwilowo
1443 // Error err; // MD
1444 // err<<"Jeszcze nie ma normay Airy'ego\n";
1445 // throw err;
1446 wartosc = airy_kwadr_pierwotna(x_kon, E, A, B) - airy_kwadr_pierwotna(x_pocz, E, A, B);
1447 }
1448 else
1449 {
1450 if(E >= y_pocz)
1451 {
1452 wartosc = tryg_kwadr_pierwotna(x_kon, E, A, B) - tryg_kwadr_pierwotna(x_pocz, E, A, B);
1453 }
1454 else
1455 {
1456 wartosc = exp_kwadr_pierwotna(x_kon, E, A, B) - exp_kwadr_pierwotna(x_pocz, E, A, B);
1457 }
1458 }
1459 return wartosc;
1460}
1461/*****************************************************************************/
1462 double warstwa::Eodk(double k) const { return k * k / (2 * masa_r); }
1463/*****************************************************************************/
1464 void warstwa_skraj::przesun_igreki(double dE)
1465{
1466 y += dE;
1467 warstwa::przesun_igreki(dE);
1468}
1469/*****************************************************************************/
1470 double warstwa_skraj::ffala(double x, double E) const
1471{
1472 double wartosc;
1473 if(lp == lewa)
1474 {
1475 wartosc = 0; // Po lewej nie ma funkcji a
1476 }
1477 else
1478 {
1479 if(E > y)
1480 {
1481 Error err; // MD
1482 err << "Energia nad skrajna bariera!\nE = " << E << " y = " << y << "\n";
1483 throw err;
1484 }
1485 else
1486 {
1487 // std::clog<<"Lewa, expb z x = "<<x<<" a x_warstwy = "<<(this->iks)<<"\n";
1488 wartosc = expa(x, E);
1489 }
1490 }
1491 return wartosc;
1492}
1493/*****************************************************************************/
1494 double warstwa_skraj::ffalb(double x, double E) const
1495{
1496 double wartosc;
1497 if(lp == prawa)
1498 {
1499 wartosc = 0; // Po prawej nie ma funkcji b
1500 }
1501 else
1502 {
1503 if(E > y)
1504 {
1505 Error err; // MD
1506 err << "Energia nad skrajna bariera!\nE = " << E << " y = " << y << "\n";
1507 throw err;
1508 }
1509 else
1510 {
1511 wartosc = expb(x, E);
1512 }
1513 }
1514 return wartosc;
1515}
1516/*****************************************************************************/
1517 double warstwa_skraj::ffala_prim(double x, double E) const
1518{
1519 double wartosc;
1520 if(lp == lewa)
1521 {
1522 wartosc = 0; // Po lewej nie ma funkcji a
1523 }
1524 else
1525 {
1526 if(E > y)
1527 {
1528 Error err; // MD
1529 err << "Energia nad skrajna bariera!\nE = " << E << " y = " << y << "\n";
1530 throw err;
1531 }
1532 else
1533 {
1534 wartosc = expa_prim(x, E);
1535 }
1536 }
1537 return wartosc;
1538}
1539/*****************************************************************************/
1540 double warstwa_skraj::ffalb_prim(double x, double E) const
1541{
1542 double wartosc;
1543 if(lp == prawa)
1544 {
1545 wartosc = 0; // Po lewej nie ma funkcji b
1546 }
1547 else
1548 {
1549 if(E > y)
1550 {
1551 Error err; // MD
1552 err << "Energia nad skrajna bariera!\n";
1553 throw err;
1554 }
1555 else
1556 {
1557 wartosc = expb_prim(x, E);
1558 }
1559 }
1560 return wartosc;
1561}
1562/*****************************************************************************/
1563 int warstwa_skraj::zera_ffal(double, double, double) const { return 0; }
1564/*****************************************************************************/
1565 double warstwa_skraj::norma_kwadr(double E, double C) const
1566{
1567 if(E > y)
1568 {
1569 Error err; // MD
1570 err << "Zla energia!\n";
1571 throw err;
1572 }
1573 double kp = sqrt(2 * masa_p * (y - E));
1574 return C * C / (2 * kp);
1575}
1576/*****************************************************************************/
1577 double warstwa_skraj::funkcjafal(double x, double E, double C) const
1578{
1579 double wynik;
1580 if(lp == lewa)
1581 {
1582 wynik = C * ffalb(x, E);
1583 }
1584 else
1585 {
1586 wynik = C * ffala(x, E);
1587 }
1588 return wynik;
1589}
1590/*****************************************************************************/
1591 double warstwa_skraj::funkcjafal_prim(double x, double E, double C) const
1592{
1593 double wynik;
1594 if(lp == lewa)
1595 {
1596 wynik = C * ffalb_prim(x, E);
1597 }
1598 else
1599 {
1600 wynik = C * ffala_prim(x, E);
1601 }
1602 return wynik;
1603}
1604/*****************************************************************************/
1605
1606/*****************************************************************************/
1607 punkt::punkt() {}
1608/*****************************************************************************/
1609 punkt::punkt(double e, double w)
1610{
1611 en = e;
1612 wart = w;
1613}
1614/*****************************************************************************/
1615 punkt::punkt(const stan & st)
1616{
1617 en = st.poziom;
1618 wart = 0;
1619}
1620/*****************************************************************************/
1621
1622/*****************************************************************************/
1623 stan::stan() { liczba_zer = -1; }
1624/*****************************************************************************/
1625 stan::stan(double E, A2D & V, int lz)
1626{
1627 poziom = E;
1628 int M = V.dim1();
1629 wspolczynniki.resize(M);
1630 for(int i = 0; i <= M - 1; i++)
1631 {
1632 wspolczynniki[i] = V[i][M - 1];
1633 }
1634 liczba_zer = lz;
1635 prawdopodobienstwa.reserve(M / 2 + 1);
1636}
1637/*****************************************************************************/
1638 void stan::przesun_poziom(double dE) { poziom += dE; }
1639/*****************************************************************************/
1640
1641/***************************************************************************** Stara wersja
1642 struktura::struktura(const std::vector<warstwa> & tablica)
1643{
1644 gora = tablica[0].y_kon;
1645 dol = gora;
1646 double czydol;
1647 if( tablica[0].pole != 0 || tablica[tablica.size() - 1].pole != 0 || tablica[0].y_kon != tablica[tablica.size() -
1648 1].y_pocz)
1649 {
1650 Error err; // MD
1651 err<<"Zle energie skajnych warstw!\n";
1652 throw err;
1653 }
1654 kawalki.push_back(tablica.front());
1655 for(int i = 1; i <= (int) tablica.size() - 2; i++)
1656 {
1657 kawalki.push_back(tablica[i]);
1658 czydol = (tablica[i].y_pocz > tablica[i].y_kon)?tablica[i].y_kon:tablica[i].y_pocz;
1659 if(czydol < dol)
1660 {
1661 dol = czydol;
1662 }
1663 if(tablica[i].pole == 0)
1664 {
1665 progi.push_back(tablica[i].y_pocz);
1666 }
1667 }
1668 kawalki.push_back(tablica.back());
1669 if(dol >= gora)
1670 {
1671 Error err; // MD
1672 err<<"Brak jakiejkolwiek studni!\n";
1673 throw err;
1674 }
1675 std::vector<double>::iterator it = progi.begin();
1676 while(it != progi.end())
1677 {
1678 if( *it == dol)
1679 {
1680 progi.erase(it);
1681 }
1682 it++;
1683 }
1684 dokl = 1e-6;
1685}
1686*****************************************************************************/
1687 struktura::struktura(const std::vector<warstwa *> & tablica, rodzaj co)
1688{
1689 lewa = *((const warstwa_skraj *)tablica[0]);
1690 if(lewa.lp == warstwa_skraj::prawa)
1691 {
1692 Error err; // MD
1693 err << "Pierwsza warstwa nie jest lewa!\n";
1694 throw err;
1695 }
1696 gora = lewa.y; // Zero to warstwa skrajna
1697 dol = gora;
1698 prawa = *((const warstwa_skraj *)tablica.back());
1699 if(prawa.lp == warstwa_skraj::lewa)
1700 {
1701 Error err; // MD
1702 err << "Ostatnia warstwa nie jest prawa!\n";
1703 throw err;
1704 }
1705
1706 double czydol;
1707 if(lewa.y != prawa.y)
1708 {
1709 Error err; // MD
1710 err << "Zle energie skajnych warstw!\n";
1711 throw err;
1712 }
1713 int i;
1714 for(i = 1; i <= (int)tablica.size() - 2; i++)
1715 {
1716 if(tablica[i - 1]->x_kon != tablica[i]->x_pocz)
1717 {
1718 Error err; // MD
1719 err << "Rozne krance warstw " << (i - 1) << " i " << i << " w " << co << ": " << (tablica[i - 1]->x_kon)
1720 << " =/= " << (tablica[i]->x_pocz) << "\n";
1721 throw err;
1722 }
1723 kawalki.push_back(*tablica[i]);
1724 tablica[i - 1]->nast = tablica[i]; // ustawianie wskaznika na sasiadke
1725 czydol = (tablica[i]->y_pocz > tablica[i]->y_kon) ? tablica[i]->y_kon : tablica[i]->y_pocz;
1726 if(czydol < dol)
1727 {
1728 dol = czydol;
1729 }
1730 if(tablica[i]->pole == 0)
1731 {
1732 progi.push_back(tablica[i]->y_pocz);
1733 }
1734 }
1735 // std::clog<<"i = "<<i<<"\tx_pocz("<<(i-1)"<<(tablica[i - 1]->x_pocz)<<"\n";
1736 if(tablica[i - 1]->x_kon != tablica[i]->x_pocz)
1737 {
1738 Error err; // MD
1739 err << "Rozne krance warstw " << (i - 1) << " i " << i << "\n";
1740 throw err;
1741 }
1742 if(dol >= gora)
1743 {
1744 Error err; // MD
1745 err << "Brak jakiejkolwiek studni!\n";
1746 throw err;
1747 }
1748 double gornyprog = dol; // pierwsze prog _pod_ bariera
1749 std::vector<double>::iterator it = progi.begin();
1750 while(it != progi.end())
1751 {
1752 // std::clog<<"prog = "<<(*it)<<"\n";
1753 if(*it > gornyprog && *it < 0.)
1754 {
1755 gornyprog = *it;
1756 }
1757 if(*it == dol)
1758 {
1759 it = progi.erase(it);
1760 }
1761 else
1762 {
1763 it++;
1764 }
1765 }
1766 typ = co;
1767 dokl = 1e-6;
1768 double start;
1769 start = -dokl; //(0.001*gornyprog > -dokl)?(0.001*gornyprog):-dokl; // zeby nie zaczac od progu
1770 szukanie_poziomow(start);
1771 normowanie();
1772 // profil(0., 1e-5);
1773}
1774/*****************************************************************************/
1775#ifdef MICHAL
1776 struktura::struktura(std::ifstream & plik, rodzaj co, bool bezlicz)
1777{
1778 std::string wiersz, bezkoment;
1779 boost::regex wykoment("#.*");
1780 std::string nic("");
1781 boost::regex pust("\\s+");
1782 boost::regex pole("pole");
1783 std::vector<double> parametry;
1784 double liczba;
1785 double x_pocz = 0, x_kon, y_pocz, y_kon, npar1, npar2, masa_p, masa_r;
1786 bool bylalewa = false, bylawew = false, bylaprawa = false;
1787 std::vector<warstwa *> tablica;
1788 warstwa * wskazwar;
1789 std::getline(plik, wiersz);
1790 bool jestpole = regex_search(wiersz, pole);
1791 bool gwiazdka = false;
1792 std::clog << "\njestpole = " << jestpole << "\n";
1793 int max_par = (jestpole) ? 7 : 6;
1794 int min_par = (jestpole) ? 5 : 4; // maksymalna i minimalna liczba parametrow w wierszu
1795 while(!plik.eof())
1796 {
1797 bezkoment = regex_replace(wiersz, wykoment, nic);
1798 boost::sregex_token_iterator it(bezkoment.begin(), bezkoment.end(), pust, -1);
1799 boost::sregex_token_iterator kon; // z jakiegos powodu oznacza to koniec
1800 if(it != kon) // niepusta zawartosc
1801 {
1802 parametry.clear();
1803 while(it != kon)
1804 {
1805 std::clog << *it << " * ";
1806 if(it->str() == "*") // oznaczenie warstw do liczenia nnosnikow
1807 {
1808 gwiazdka = true;
1809 ++it;
1810 continue;
1811 }
1812 try
1813 {
1814 liczba = boost::lexical_cast<double>(*it);
1815 it++;
1816 }
1817 catch(boost::bad_lexical_cast &)
1818 {
1819 Error err; // MD
1820 err << "\n napis " << *it << " nie jest liczba\n";
1821 throw err;
1822 }
1823 parametry.push_back(liczba);
1824 std::clog << "\n";
1825 }
1826 if(bylalewa &&
1827 (((int)parametry.size() < min_par && parametry.size() != 1) || (int)parametry.size() > max_par))
1828 {
1829 Error err; // MD
1830 if(jestpole)
1831 err << "\nwarstwa wymaga 1 lub od 5 do 7 parametrow, a sa " << parametry.size() << "\n";
1832 else
1833 err << "\nwarstwa wymaga 1 lub od 4 do 6 parametrow, a sa " << parametry.size() << "\n";
1834 throw err;
1835 }
1836 if(bylalewa)
1837 {
1838 if(parametry.size() == 1) // prawa
1839 {
1840 bylaprawa = true;
1841 wskazwar = new warstwa_skraj(warstwa_skraj::prawa, parametry[0], parametry[0], x_pocz, 0.);
1842 std::clog << "\nrobi sie prawa: masa = " << parametry[0] << " x_pocz = " << x_pocz << "\n";
1843 tablica.push_back(wskazwar);
1844 break;
1845 }
1846 else // wewnetrzna
1847 {
1848 x_kon = x_pocz + parametry[0];
1849 y_pocz = -parametry[1];
1850 if(jestpole)
1851 {
1852 y_kon = -parametry[2];
1853 ;
1854 masa_p = parametry[3];
1855 masa_r = parametry[4];
1856 }
1857 else
1858 {
1859 y_kon = y_pocz;
1860 masa_p = parametry[2];
1861 masa_r = parametry[3];
1862 }
1863 npar1 = 0.;
1864 npar2 = 0.;
1865 if((parametry[0] <= 0.) || (masa_p <= 0.) || (masa_r <= 0.))
1866 {
1867 Error err; // MD
1868 err << "\nAaaaaa! x_kon = " << parametry[0] << "masa_p = " << masa_p << "masa_r = " << masa_r
1869 << "\n";
1870 throw err;
1871 }
1872 bylalewa = true;
1873 if((int)parametry.size() == min_par + 1)
1874 {
1875 npar1 = parametry[min_par];
1876 }
1877 if((int)parametry.size() == min_par + 2)
1878 {
1879 npar1 = parametry[min_par];
1880 npar2 = parametry[min_par + 1];
1881 }
1882 wskazwar = new warstwa(masa_p, masa_r, x_pocz, y_pocz, x_kon, y_kon, npar1, npar2);
1883 std::clog << "masa_p = " << masa_p << ", masa_r = " << masa_r << ", x_pocz = " << x_pocz
1884 << ", y_pocz = " << y_pocz << ", x_kon = " << x_kon << ", y_kon = " << y_kon
1885 << ", npar1 = " << npar1 << ", npar2 = " << npar2 << "\n";
1886 ;
1887 tablica.push_back(wskazwar);
1888 if(gwiazdka)
1889 {
1890 gwiazdki.push_back(tablica.size() - 2);
1891 gwiazdka = false;
1892 }
1893 x_pocz = x_kon;
1894 }
1895 }
1896 else
1897 {
1898 bylalewa = true;
1899 wskazwar = new warstwa_skraj(warstwa_skraj::lewa, parametry[0], parametry[0], x_pocz, 0.);
1900 tablica.push_back(wskazwar);
1901 std::clog << "\nrobi sie lewa: masa = " << parametry[0] << " x_pocz = " << x_pocz << "\n";
1902 }
1903 }
1904 if(bylalewa && bylawew && bylaprawa)
1905 {
1906 std::clog << "\nWsystko bylo\n";
1907 }
1908 std::getline(plik, wiersz);
1909 }
1910
1911 // ponizej zawartosc konstruktora od tablicy
1912
1913 lewa = *((const warstwa_skraj *)tablica[0]);
1914 if(lewa.lp == warstwa_skraj::prawa)
1915 {
1916 Error err; // MD
1917 err << "Pierwsza warstwa nie jest lewa!\n";
1918 throw err;
1919 }
1920 gora = lewa.y; // Zero to warstwa skrajna
1921 dol = gora;
1922 prawa = *((const warstwa_skraj *)tablica.back());
1923 if(prawa.lp == warstwa_skraj::lewa)
1924 {
1925 Error err; // MD
1926 err << "Ostatnia warstwa nie jest prawa!\n";
1927 throw err;
1928 }
1929
1930 double czydol;
1931 if(lewa.y != prawa.y)
1932 {
1933 Error err; // MD
1934 err << "Zle energie skajnych warstw!\n";
1935 throw err;
1936 }
1937 int i;
1938 for(i = 1; i <= (int)tablica.size() - 2; i++)
1939 {
1940 if(tablica[i - 1]->x_kon != tablica[i]->x_pocz)
1941 {
1942 Error err; // MD
1943 err << "Rozne krance warstw " << (i - 1) << " i " << i << " w " << co << ": " << (tablica[i - 1]->x_kon)
1944 << " =/= " << (tablica[i]->x_pocz) << "\n";
1945 throw err;
1946 }
1947 kawalki.push_back(*tablica[i]);
1948 tablica[i - 1]->nast = tablica[i]; // ustawianie wskaznika na sasiadke
1949 czydol = (tablica[i]->y_pocz > tablica[i]->y_kon) ? tablica[i]->y_kon : tablica[i]->y_pocz;
1950 if(czydol < dol)
1951 {
1952 dol = czydol;
1953 }
1954 if(tablica[i]->pole == 0)
1955 {
1956 progi.push_back(tablica[i]->y_pocz);
1957 }
1958 }
1959 // std::clog<<"i = "<<i<<"\tx_pocz("<<(i-1)"<<(tablica[i - 1]->x_pocz)<<"\n";
1960 if(tablica[i - 1]->x_kon != tablica[i]->x_pocz)
1961 {
1962 Error err; // MD
1963 err << "Rozne krance warstw " << (i - 1) << " i " << i << "\n";
1964 throw err;
1965 }
1966 if(dol >= gora)
1967 {
1968 Error err; // MD
1969 err << "Brak jakiejkolwiek studni!\n";
1970 throw err;
1971 }
1972 double gornyprog = dol; // pierwsze prog _pod_ bariera
1973 // std::clog<<"gornyprog = "<<gornyprog<<"\n";
1974 std::vector<double>::iterator it = progi.begin();
1975 while(it != progi.end())
1976 {
1977 // std::clog<<"prog = "<<(*it)<<"\n";
1978 if(*it > gornyprog && *it < 0.)
1979 {
1980 gornyprog = *it;
1981 }
1982 if(*it == dol)
1983 {
1984 it = progi.erase(it);
1985 }
1986 else
1987 {
1988 it++;
1989 }
1990 }
1991 typ = co;
1992 if(co == el)
1993 nazwa = "el";
1994 if(co == hh)
1995 nazwa = "hh";
1996 if(co == lh)
1997 nazwa = "lh";
1998
1999 if(!bezlicz)
2000 {
2001 //1.5
2002 // dokl = 1e-7;
2003 dokl = 1e-10;
2004 //
2005 double start;
2006 // std::clog<<"gornyprog = "<<gornyprog<<"\n";
2007 start = -dokl; //(0.001*gornyprog > -dokl)?(0.001*gornyprog):-dokl; // zeby nie zaczac od progu
2008 // szukanie_poziomow_2(start,dokl,false);
2009 // std::cout<<"\nTeraz z zpoprawka2\n"<<std::endl;
2010 szukanie_poziomow_zpoprawka2(start,dokl);
2011 normowanie();
2012 }
2013}
2014/*****************************************************************************/
2015 struktura::struktura(std::ifstream & plik, const std::vector<double> & poziomy, rodzaj co)
2016{
2017 std::string wiersz, bezkoment;
2018 boost::regex wykoment("#.*");
2019 std::string nic("");
2020 boost::regex pust("\\s+");
2021 boost::regex pole("pole");
2022 std::vector<double> parametry;
2023 double liczba;
2024 double x_pocz = 0, x_kon, y_pocz, y_kon, npar1, npar2, masa_p, masa_r;
2025 bool bylalewa = false, bylawew = false, bylaprawa = false;
2026 std::vector<warstwa *> tablica;
2027 warstwa * wskazwar;
2028 std::getline(plik, wiersz);
2029 bool jestpole = regex_search(wiersz, pole);
2030 bool gwiazdka = false;
2031 std::clog << "\njestpole = " << jestpole << "\n";
2032 int max_par = (jestpole) ? 7 : 6;
2033 int min_par = (jestpole) ? 5 : 4; // maksymalna i minimalna liczba parametrow w wierszu
2034 while(!plik.eof())
2035 {
2036 bezkoment = regex_replace(wiersz, wykoment, nic);
2037 boost::sregex_token_iterator it(bezkoment.begin(), bezkoment.end(), pust, -1);
2038 boost::sregex_token_iterator kon; // z jakiegos powodu oznacza to koniec
2039 if(it != kon) // niepusta zawartosc
2040 {
2041 parametry.clear();
2042 while(it != kon)
2043 {
2044 std::clog << *it << " * ";
2045 if(it->str() == "*") // oznaczenie warstw do liczenia nnosnikow
2046 {
2047 gwiazdka = true;
2048 ++it;
2049 continue;
2050 }
2051 try
2052 {
2053 liczba = boost::lexical_cast<double>(*it);
2054 it++;
2055 }
2056 catch(boost::bad_lexical_cast &)
2057 {
2058 Error err; // MD
2059 err << "\n napis " << *it << " nie jest liczba\n";
2060 throw err;
2061 }
2062 parametry.push_back(liczba);
2063 std::clog << "\n";
2064 }
2065 if(bylalewa &&
2066 (((int)parametry.size() < min_par && parametry.size() != 1) || (int)parametry.size() > max_par))
2067 {
2068 Error err; // MD
2069 if(jestpole)
2070 std::cerr << "\nwarstwa wymaga 1 lub od 5 do 7 parametrow, a sa " << parametry.size() << "\n";
2071 else
2072 std::cerr << "\nwarstwa wymaga 1 lub od 4 do 6 parametrow, a sa " << parametry.size() << "\n";
2073 throw err;
2074 }
2075 if(bylalewa)
2076 {
2077 if(parametry.size() == 1) // prawa
2078 {
2079 bylaprawa = true;
2080 wskazwar = new warstwa_skraj(warstwa_skraj::prawa, parametry[0], parametry[0], x_pocz, 0.);
2081 std::clog << "\nrobi sie prawa: masa = " << parametry[0] << " x_pocz = " << x_pocz << "\n";
2082 tablica.push_back(wskazwar);
2083 break;
2084 }
2085 else // wewnetrzna
2086 {
2087 x_kon = x_pocz + parametry[0];
2088 y_pocz = -parametry[1];
2089 if(jestpole)
2090 {
2091 y_kon = -parametry[2];
2092 ;
2093 masa_p = parametry[3];
2094 masa_r = parametry[4];
2095 }
2096 else
2097 {
2098 y_kon = y_pocz;
2099 masa_p = parametry[2];
2100 masa_r = parametry[3];
2101 }
2102 npar1 = 0.;
2103 npar2 = 0.;
2104 if((parametry[0] <= 0.) || (masa_p <= 0.) || (masa_r <= 0.))
2105 {
2106 Error err; // MD
2107 err << "\nAaaaaa! x_kon = " << parametry[0] << "masa_p = " << masa_p << "masa_r = " << masa_r
2108 << "\n";
2109 throw err;
2110 }
2111 bylalewa = true;
2112 if((int)parametry.size() == min_par + 1)
2113 {
2114 npar1 = parametry[min_par];
2115 }
2116 if((int)parametry.size() == min_par + 2)
2117 {
2118 npar1 = parametry[min_par];
2119 npar2 = parametry[min_par + 1];
2120 }
2121 wskazwar = new warstwa(masa_p, masa_r, x_pocz, y_pocz, x_kon, y_kon, npar1, npar2);
2122 std::clog << "masa_p = " << masa_p << ", masa_r = " << masa_r << ", x_pocz = " << x_pocz
2123 << ", y_pocz = " << y_pocz << ", x_kon = " << x_kon << ", y_kon = " << y_kon
2124 << ", npar1 = " << npar1 << ", npar2 = " << npar2 << "\n";
2125 ;
2126 tablica.push_back(wskazwar);
2127 if(gwiazdka)
2128 {
2129 gwiazdki.push_back(tablica.size() - 2);
2130 gwiazdka = false;
2131 }
2132 x_pocz = x_kon;
2133 }
2134 }
2135 else
2136 {
2137 bylalewa = true;
2138 wskazwar = new warstwa_skraj(warstwa_skraj::lewa, parametry[0], parametry[0], x_pocz, 0.);
2139 tablica.push_back(wskazwar);
2140 std::clog << "\nrobi sie lewa: masa = " << parametry[0] << " x_pocz = " << x_pocz << "\n";
2141 }
2142 }
2143 if(bylalewa && bylawew && bylaprawa)
2144 {
2145 std::clog << "\nWsystko bylo\n";
2146 }
2147 std::getline(plik, wiersz);
2148 }
2149
2150 // ponizej zawartosc konstruktora od tablicy
2151
2152 lewa = *((const warstwa_skraj *)tablica[0]);
2153 if(lewa.lp == warstwa_skraj::prawa)
2154 {
2155 Error err; // MD
2156 err << "Pierwsza warstwa nie jest lewa!\n";
2157 throw err;
2158 }
2159 gora = lewa.y; // Zero to warstwa skrajna
2160 dol = gora;
2161 prawa = *((const warstwa_skraj *)tablica.back());
2162 if(prawa.lp == warstwa_skraj::lewa)
2163 {
2164 Error err; // MD
2165 err << "Ostatnia warstwa nie jest prawa!\n";
2166 throw err;
2167 }
2168
2169 double czydol;
2170 if(lewa.y != prawa.y)
2171 {
2172 Error err; // MD
2173 err << "Zle energie skajnych warstw!\n";
2174 throw err;
2175 }
2176 int i;
2177 for(i = 1; i <= (int)tablica.size() - 2; i++)
2178 {
2179 if(tablica[i - 1]->x_kon != tablica[i]->x_pocz)
2180 {
2181 Error err; // MD
2182 err << "Rozne krance warstw " << (i - 1) << " i " << i << " w " << co << ": " << (tablica[i - 1]->x_kon)
2183 << " =/= " << (tablica[i]->x_pocz) << "\n";
2184 throw err;
2185 }
2186 kawalki.push_back(*tablica[i]);
2187 tablica[i - 1]->nast = tablica[i]; // ustawianie wskaznika na sasiadke
2188 czydol = (tablica[i]->y_pocz > tablica[i]->y_kon) ? tablica[i]->y_kon : tablica[i]->y_pocz;
2189 if(czydol < dol)
2190 {
2191 dol = czydol;
2192 }
2193 if(tablica[i]->pole == 0)
2194 {
2195 progi.push_back(tablica[i]->y_pocz);
2196 }
2197 }
2198 // std::clog<<"i = "<<i<<"\tx_pocz("<<(i-1)"<<(tablica[i - 1]->x_pocz)<<"\n";
2199 if(tablica[i - 1]->x_kon != tablica[i]->x_pocz)
2200 {
2201 Error err; // MD
2202 err << "Rozne krance warstw " << (i - 1) << " i " << i << "\n";
2203 throw err;
2204 }
2205 if(dol >= gora)
2206 {
2207 Error err; // MD
2208 err << "Brak jakiejkolwiek studni!\n";
2209 throw err;
2210 }
2211 double gornyprog = dol; // pierwsze prog _pod_ bariera
2212 // std::clog<<"gornyprog = "<<gornyprog<<"\n";
2213 std::vector<double>::iterator it = progi.begin();
2214 while(it != progi.end())
2215 {
2216 // std::clog<<"prog = "<<(*it)<<"\n";
2217 if(*it > gornyprog && *it < 0.)
2218 {
2219 gornyprog = *it;
2220 }
2221 if(*it == dol)
2222 {
2223 it = progi.erase(it);
2224 }
2225 else
2226 {
2227 it++;
2228 }
2229 }
2230 typ = co;
2231 if(co == el)
2232 nazwa = "el";
2233 if(co == hh)
2234 nazwa = "hh";
2235 if(co == lh)
2236 nazwa = "lh";
2237
2238 stany_z_tablicy(poziomy);
2239 normowanie();
2240}
2241#endif
2242/*****************************************************************************/
2243 void struktura::przesun_energie(double dE)
2244{
2245 gora += dE;
2246 dol += dE;
2247 lewa.przesun_igreki(dE);
2248 prawa.przesun_igreki(dE);
2249 int i; // LUKASZ Visual 6
2250 for(/*int*/ i = 0; i <= (int)kawalki.size() - 1; i++) // LUKASZ Visual 6
2251 {
2252 kawalki[i].przesun_igreki(dE);
2253 }
2254 for(/*int*/ i = 0; i <= (int)progi.size() - 1; i++) // LUKASZ Visual 6
2255 {
2256 progi[i] += dE;
2257 }
2258 for(/*int*/ i = 0; i <= (int)rozwiazania.size() - 1; i++) // LUKASZ Visual 6
2259 {
2260 rozwiazania[i].przesun_poziom(dE);
2261 }
2262}
2263/*****************************************************************************/
2264 double struktura::dlugosc_z_A(const double dlugA) { return dlugA / przelm; }
2265/*****************************************************************************/
2266 double struktura::dlugosc_na_A(const double dlug) { return dlug * przelm; }
2267/*****************************************************************************/
2268 double struktura::koncentracja_na_cm_3(const double k_w_wew) { return k_w_wew / (przelm * przelm * przelm) * 1e24; }
2269/*****************************************************************************/
2270 double struktura::czyosobliwa(double E)
2271{
2272 int N = kawalki.size() + 2; // liczba warstw
2273 // Bylo bez '+ 2'
2274 if(N < 3)
2275 {
2276 Error err; // MD
2277 err << "Za malo warstw, bo " << N << "\n";
2278 throw err;
2279 }
2280 int M = 2 * N - 2; // liczba rownan
2281 A2D macierz(M, M, 0.0);
2282 zrobmacierz(E, macierz);
2283 // std::clog<<E<<"\n"<<macierz<<"\n";
2284 A1D S(macierz.dim1());
2285 JAMA::SVD<double> rozklad(macierz);
2286 rozklad.getSingularValues(S);
2287
2288 A2D V(M, M);
2289 A2D U(M, M);
2290 rozklad.getV(V);
2291 rozklad.getU(U);
2292 A2D UV = matmult(U, V);
2293 // std::clog<<E<<"\nV = "<<V<<"\n";
2294 JAMA::LU<double> rozkladUV(UV);
2295 double detUV = rozkladUV.det();
2296 // return((S[S.dim() - 1]));
2297
2298 // return((S[S.dim() - 1])/S[S.dim() - 2]); // Dzielone, zeby nie bylo zer podwojnych
2299 double dzielnik = 1; // Zeby wyeleminowac zera na plaskich kawalkach
2300 for(int i = 0; i <= (int)progi.size() - 1; i++)
2301 {
2302 if(E == progi[i])
2303 {
2304 Error err; // MD
2305 err << "Energia " << E << " rowna progowi nr " << i << "\n";
2306 throw err;
2307 }
2308 dzielnik *= (E - progi[i]);
2309 }
2310 return (detUV * S[S.dim() - 1]) / dzielnik;
2311 // return(S[S.dim() - 1]);
2312}
2313/***************************************************************************** Stare
2314 void struktura::zrobmacierz(double E, A2D & macierz)
2315{
2316 int N = kawalki.size(); // liczba warstw
2317 double x = kawalki[1].x_pocz;
2318 macierz[0][0] = kawalki[0].ffalb(x, E);
2319 macierz[0][1] = -kawalki[1].ffala(x, E);
2320 macierz[0][2] = -kawalki[1].ffalb(x, E);
2321 macierz[1][0] = kawalki[0].ffalb_prim(x, E)/kawalki[0].masa;
2322 macierz[1][1] = -kawalki[1].ffala_prim(x, E)/kawalki[1].masa;
2323 macierz[1][2] = -kawalki[1].ffalb_prim(x, E)/kawalki[1].masa;
2324 int n = 1;
2325 for(n = 1; n <= N - 3; n++)
2326 {
2327 x = kawalki[n + 1].x_pocz;
2328 macierz[2*n][2*n - 1] = kawalki[n].ffala(x, E);
2329 macierz[2*n][2*n] = kawalki[n].ffalb(x, E);
2330 macierz[2*n][2*n + 1] = -kawalki[n+1].ffala(x, E);
2331 macierz[2*n][2*n + 2] = -kawalki[n+1].ffalb(x, E);
2332
2333 macierz[2*n + 1][2*n - 1] = kawalki[n].ffala_prim(x, E)/kawalki[n].masa;
2334 macierz[2*n + 1][2*n] = kawalki[n].ffalb_prim(x, E)/kawalki[n].masa;
2335 macierz[2*n + 1][2*n + 1] = -kawalki[n+1].ffala_prim(x, E)/kawalki[n + 1].masa;
2336 macierz[2*n + 1][2*n + 2] = -kawalki[n+1].ffalb_prim(x, E)/kawalki[n + 1].masa;
2337 }
2338 x = kawalki[N - 1].x_pocz;
2339 // std::clog<"ostatni x = "<<(x*warstwa::przelm)<<"\n";
2340 macierz[2*n][2*n - 1] = kawalki[n].ffala(x, E);
2341 macierz[2*n][2*n] = kawalki[n].ffalb(x, E);
2342 macierz[2*n][2*n + 1] = -kawalki[n+1].ffala(x, E);
2343
2344 macierz[2*n + 1][2*n - 1] = kawalki[n].ffala_prim(x, E)/kawalki[n].masa;
2345 macierz[2*n + 1][2*n] = kawalki[n].ffalb_prim(x, E)/kawalki[n].masa;
2346 macierz[2*n + 1][2*n + 1] = -kawalki[n+1].ffala_prim(x, E)/kawalki[n + 1].masa;
2347 }
2348*****************************************************************************/
2349 void struktura::zrobmacierz(double E, A2D & macierz)
2350{
2351 int N = kawalki.size() + 2; // liczba warstw
2352 double x = lewa.iks;
2353 macierz[0][0] = lewa.ffalb(x, E);
2354 macierz[0][1] = -kawalki[0].ffala(x, E);
2355 macierz[0][2] = -kawalki[0].ffalb(x, E);
2356 macierz[1][0] = lewa.ffalb_prim(x, E) / lewa.masa_p;
2357 macierz[1][1] = -kawalki[0].ffala_prim(x, E) / kawalki[0].masa_p(E);
2358 macierz[1][2] = -kawalki[0].ffalb_prim(x, E) / kawalki[0].masa_p(E);
2359 int n = 1;
2360 for(n = 1; n <= N - 3; n++)
2361 {
2362 x = kawalki[n].x_pocz;
2363 macierz[2 * n][2 * n - 1] = kawalki[n - 1].ffala(x, E);
2364 macierz[2 * n][2 * n] = kawalki[n - 1].ffalb(x, E);
2365 macierz[2 * n][2 * n + 1] = -kawalki[n].ffala(x, E);
2366 macierz[2 * n][2 * n + 2] = -kawalki[n].ffalb(x, E);
2367
2368 macierz[2 * n + 1][2 * n - 1] = kawalki[n - 1].ffala_prim(x, E) / kawalki[n - 1].masa_p(E);
2369 macierz[2 * n + 1][2 * n] = kawalki[n - 1].ffalb_prim(x, E) / kawalki[n - 1].masa_p(E);
2370 macierz[2 * n + 1][2 * n + 1] = -kawalki[n].ffala_prim(x, E) / kawalki[n].masa_p(E);
2371 macierz[2 * n + 1][2 * n + 2] = -kawalki[n].ffalb_prim(x, E) / kawalki[n].masa_p(E);
2372 }
2373 x = prawa.iks;
2374 // std::clog<<"ostatni x = "<<(x*warstwa::przelm)<<"\n";
2375 macierz[2 * n][2 * n - 1] = kawalki[n - 1].ffala(x, E);
2376 macierz[2 * n][2 * n] = kawalki[n - 1].ffalb(x, E);
2377 macierz[2 * n][2 * n + 1] = -prawa.ffala(x, E);
2378
2379 macierz[2 * n + 1][2 * n - 1] = kawalki[n - 1].ffala_prim(x, E) / kawalki[n - 1].masa_p(E);
2380 macierz[2 * n + 1][2 * n] = kawalki[n - 1].ffalb_prim(x, E) / kawalki[n - 1].masa_p(E);
2381 macierz[2 * n + 1][2 * n + 1] = -prawa.ffala_prim(x, E) / prawa.masa_p;
2382}
2383/*****************************************************************************/
2384 double struktura::sprawdz_ciaglosc(double E, A2D & V) // zwraca sume skokow funkcji na interfejsach. znak ujemny oznacza
2385 // ze byla zmiana znaku na nieciaNlosci na interfejsie
2386{
2387 int N = kawalki.size() + 2; // liczba warstw
2388 int M = 2 * N - 2; // liczba rownan
2389 A2D macierz(M, M, 0.0);
2390 zrobmacierz(E, macierz);
2391 JAMA::SVD<double> rozklad(macierz);
2392 rozklad.getV(V);
2393 double Al, Bl, A, B;
2394 double lewawart, prawawart;
2395 double sumaroz = 0;
2396 double znak = 1.;
2397 Al = V[0][V.dim2() - 1];
2398 A = V[1][V.dim2() - 1];
2399 B = V[2][V.dim2() - 1];
2400 lewawart = lewa.funkcjafal(lewa.iks, E, Al);
2401 prawawart = kawalki[0].funkcjafal(kawalki[0].x_pocz, E, A, B);
2402 if(lewawart * prawawart < 0)
2403 {
2404 znak = -1.;
2405#ifdef LOGUJ // MD
2406 std::clog << "\nE = " << E << " zmiana znaku! " << lewawart << " " << prawawart << "\n";
2407#endif
2408 }
2409 sumaroz += abs(lewawart - prawawart);
2410 int j; // LUKASZ Visual 6
2411 for(/*int*/ j = 1; j <= N - 3; j++) // LUKASZ Visual 6
2412 {
2413 // std::cerr<<"j = "<<j<<" ";
2414 Al = V[2 * (j - 1) + 1][V.dim2() - 1];
2415 Bl = V[2 * (j - 1) + 2][V.dim2() - 1];
2416 A = V[2 * j + 1][V.dim2() - 1];
2417 B = V[2 * j + 2][V.dim2() - 1];
2418 lewawart = kawalki[j - 1].funkcjafal(kawalki[j - 1].x_kon, E, Al, Bl);
2419 prawawart = kawalki[j].funkcjafal(kawalki[j].x_pocz, E, A, B);
2420 if(lewawart * prawawart < 0)
2421 {
2422 znak = -1.;
2423#ifdef LOGUJ // MD
2424 std::clog << "\nE = " << E << " zmiana znaku! " << lewawart << " " << prawawart << "\n";
2425#endif
2426 }
2427 sumaroz += abs(lewawart - prawawart);
2428 }
2429 /*int*/ j = N - 2; // LUKASZ Visual 6
2430 Al = V[2 * (j - 1) + 1][V.dim2() - 1];
2431 Bl = V[2 * (j - 1) + 2][V.dim2() - 1];
2432 A = V[2 * j + 1][V.dim2() - 1];
2433 lewawart = kawalki[j - 1].funkcjafal(kawalki[j - 1].x_kon, E, Al, Bl);
2434 prawawart = prawa.funkcjafal(prawa.iks, E, A);
2435 if(lewawart * prawawart < 0)
2436 {
2437 znak = -1.;
2438#ifdef LOGUJ // MD
2439 std::clog << "\nE = " << E << " zmiana znaku! " << lewawart << " " << prawawart << "\n";
2440#endif
2441 }
2442
2443 sumaroz += abs(lewawart - prawawart);
2444 return znak * sumaroz;
2445}
2446/*****************************************************************************/
2447std::vector<std::vector<double> > struktura::rysowanie_funkcji(double E, double x0A, double xkA, double krokA)
2448{
2449 double x0 = x0A / przelm;
2450 double xk = xkA / przelm;
2451 const double krok = krokA / przelm;
2452 int N = kawalki.size() + 2; // liczba warstw
2453 int M = 2 * N - 2; // liczba rownan
2454 A2D macierz(M, M, 0.0);
2455 zrobmacierz(E, macierz);
2456 A2D V(M, M);
2457 JAMA::SVD<double> rozklad(macierz);
2458 rozklad.getV(V);
2459 // std::clog<<"To jest V:\n"<<V;
2460 double x = x0;
2461 int wsk = 0;
2462 std::vector<std::vector<double> > funkcja;
2463 funkcja.resize(2);
2464 funkcja[0].reserve(int((xk - x0) / krok));
2465 funkcja[1].reserve(int((xk - x0) / krok));
2466 double A, B;
2467 if(x < lewa.iks)
2468 {
2469 wsk = -1;
2470 }
2471 else
2472 {
2473 while((wsk <= (int)kawalki.size() - 1) && (x > kawalki[wsk].x_kon)) // Szukanie pierwszej kawalki
2474 {
2475 wsk++;
2476 }
2477 }
2478 while(x <= xk)
2479 {
2480 funkcja[0].push_back(x);
2481 if(wsk == -1.)
2482 {
2483 B = V[0][V.dim2() - 1];
2484 funkcja[1].push_back(lewa.funkcjafal(x, E, B));
2485 if(x > lewa.iks)
2486 {
2487 wsk++;
2488 }
2489 }
2490 else if(wsk <= (int)kawalki.size() - 1)
2491 {
2492 A = V[2 * wsk + 1][V.dim2() - 1];
2493 B = V[2 * wsk + 2][V.dim2() - 1];
2494 funkcja[1].push_back(kawalki[wsk].funkcjafal(x, E, A, B));
2495 }
2496 else
2497 {
2498 A = V[2 * wsk + 1][V.dim2() - 1];
2499 funkcja[1].push_back(prawa.funkcjafal(x, E, A));
2500 }
2501 x += krok;
2502 if((wsk >= 0) && (wsk <= (int)kawalki.size() - 1) && (x > kawalki[wsk].x_kon))
2503 {
2504 wsk++;
2505 }
2506 }
2507 return funkcja;
2508}
2509/*****************************************************************************
2510 int struktura::ilezer_ffal(double E)
2511{
2512 int N = kawalki.size() + 2; //liczba warstw
2513 // Bylo bez '+ 2'
2514 int M = 2*N - 2; // liczba rownan
2515 A2D macierz(M, M, 0.0);
2516 zrobmacierz(E, macierz);
2517 A2D V(M, M);
2518 JAMA::SVD<double> rozklad(macierz);
2519 rozklad.getV(V);
2520 int sumazer = 0;
2521 double A, B, As, Bs;
2522 // bool juz = false, jeszcze = true; // czy juz lub jeszcze warto sprawdzac (nie bedzie zer w warstwie, w ktorej
2523 poziom jest pod przerwa i tak samo jest we wszystkich od lub do skrajnej) int pierwsza = -1; do
2524 {
2525 pierwsza++;
2526 }while(pierwsza <= N-3 && (kawalki[pierwsza].y_pocz > E && kawalki[pierwsza].y_kon > E) );
2527 int ostatnia = N-2;
2528 do
2529 {
2530 ostatnia--;
2531 }while(ostatnia >= 0 && (kawalki[ostatnia].y_pocz > E && kawalki[ostatnia].y_kon > E) );
2532 std::clog<<"\npierwsza sprawdzana = "<<pierwsza<<" ostatnia = "<<ostatnia;
2533
2534 double sasiad; // wartosc ffal na lewym brzegu sasiada z prawej (ta, ktora powinna byc wspolna do obu)
2535 for(int j = pierwsza; j <= ostatnia - 1; j++)
2536 {
2537 A = V[2*j+1][V.dim2() - 1];
2538 B = V[2*j+2][V.dim2() - 1];
2539 As = V[2*(j+1)+1][V.dim2() - 1];
2540 Bs = V[2*(j+1)+2][V.dim2() - 1];
2541 sasiad = kawalki[j+1].funkcjafal(kawalki[j+1].x_pocz, E, As, Bs);
2542 sumazer += kawalki[j].zera_ffal(E, A, B, sasiad);
2543 }
2544 A = V[2*ostatnia+1][V.dim2() - 1];
2545 B = V[2*ostatnia+2][V.dim2() - 1];
2546 sumazer += kawalki[ostatnia].zera_ffal(E, A, B); // W ostatniej warswie nie moze byc zera na laczeniu, wiec nie ma
2547 problem return sumazer;
2548}
2549*****************************************************************************/
2550 int struktura::ilezer_ffal(double E, A2D & V)
2551{
2552 int N = kawalki.size() + 2; // liczba warstw
2553 // Bylo bez '+ 2'
2554 int M = 2 * N - 2; // liczba rownan
2555 A2D macierz(M, M, 0.0);
2556 zrobmacierz(E, macierz);
2557 JAMA::SVD<double> rozklad(macierz);
2558 rozklad.getV(V);
2559 int sumazer = 0;
2560 double A, B, Al, Bl, Ap, Bp;
2561 /*
2562 for(int i = 1; i <= N-2; i++)
2563 {
2564 A = V[2*i-1][V.dim2() - 1];
2565 B = V[2*i][V.dim2() - 1];
2566 sumazer += kawalki[i - 1].zera_ffal(E, A, B);
2567 }
2568 */ // tak bylo bez szukania pierwszej i ostatniej
2569 int pierwsza = -1;
2570 do
2571 {
2572 pierwsza++;
2573 }
2574 while(pierwsza <= N - 3 && (kawalki[pierwsza].y_pocz > E && kawalki[pierwsza].y_kon > E));
2575 int ostatnia = N - 2;
2576 do
2577 {
2578 ostatnia--;
2579 }
2580 while(ostatnia >= 0 && (kawalki[ostatnia].y_pocz > E && kawalki[ostatnia].y_kon > E));
2581 // std::clog<<"\npierwsza sprawdzana = "<<pierwsza<<" ostatnia = "<<ostatnia;
2582 double sasiadl, sasiadp; // wartosc ffal na lewym brzegu sasiada z prawej (ta, ktora powinna byc wspolna do obu)
2583 if(ostatnia == pierwsza) // tylko jedna podejrzana warstwa, nie trzeba sie sasiadami przejmowac
2584 {
2585 A = V[2 * pierwsza + 1][V.dim2() - 1];
2586 B = V[2 * pierwsza + 2][V.dim2() - 1];
2587 sumazer += kawalki[pierwsza].zera_ffal(E, A, B);
2588 }
2589 else
2590 {
2591 int j = pierwsza;
2592 A = V[2 * j + 1][V.dim2() - 1];
2593 B = V[2 * j + 2][V.dim2() - 1];
2594 Ap = V[2 * (j + 1) + 1][V.dim2() - 1];
2595 Bp = V[2 * (j + 1) + 2][V.dim2() - 1];
2596 sasiadp = kawalki[j + 1].funkcjafal(kawalki[j + 1].x_pocz, E, Ap, Bp);
2597 sasiadl = kawalki[j].funkcjafal(kawalki[j].x_pocz, E, A,
2598 B); // po lewej nie ma problemu, wiec mozna podstawic wartosc z wlasnej warstwy
2599 sumazer += kawalki[j].zera_ffal(E, A, B, sasiadl, sasiadp);
2600 for(/*int*/ j = pierwsza + 1; j <= ostatnia - 1; j++) // LUKASZ Visual 6
2601 {
2602 Al = V[2 * (j - 1) + 1][V.dim2() - 1];
2603 Bl = V[2 * (j - 1) + 2][V.dim2() - 1];
2604 A = V[2 * j + 1][V.dim2() - 1];
2605 B = V[2 * j + 2][V.dim2() - 1];
2606 Ap = V[2 * (j + 1) + 1][V.dim2() - 1];
2607 Bp = V[2 * (j + 1) + 2][V.dim2() - 1];
2608 sasiadl = kawalki[j - 1].funkcjafal(kawalki[j - 1].x_kon, E, Al, Bl);
2609 sasiadp = kawalki[j + 1].funkcjafal(kawalki[j + 1].x_pocz, E, Ap, Bp);
2610 sumazer += kawalki[j].zera_ffal(E, A, B, sasiadl, sasiadp);
2611 }
2612 j = ostatnia;
2613 A = V[2 * j + 1][V.dim2() - 1];
2614 B = V[2 * j + 2][V.dim2() - 1];
2615 Al = V[2 * (j - 1) + 1][V.dim2() - 1];
2616 Bl = V[2 * (j - 1) + 2][V.dim2() - 1];
2617 sasiadp = kawalki[j].funkcjafal(kawalki[j].x_kon, E, A,
2618 B); // po prawej nie ma problemu, wiec mozna podstawic wartosc z wlasnej warstwy
2619 sasiadl = kawalki[j - 1].funkcjafal(kawalki[j - 1].x_kon, E, Al, Bl);
2620 sumazer += kawalki[j].zera_ffal(
2621 E, A, B, sasiadl, sasiadp); // W ostatniej warswie nie moze byc zera na laczeniu, wiec nie ma problemu
2622 }
2623 return sumazer;
2624 // return 0; //do testow tylko!
2625}
2626/*****************************************************************************/
2627std::vector<double>
2628 struktura::zageszczanie(punkt p0,
2629 punkt pk) // Zageszcza az znajdzie inny znak, zaklada, ze poczatkowe znaki sa takie same
2630{
2631 std::list<punkt> lista;
2632 std::vector<double> wynik;
2633 lista.push_front(p0);
2634 lista.push_back(pk);
2635 double E;
2636 double znak = (p0.wart > 0) ? 1. : -1.;
2637 if(znak * pk.wart <= 0)
2638 {
2639 Error err; // MD
2640 err << "W zageszczaniu znaki sie roznia!\n";
2641 throw err;
2642 }
2643 std::list<punkt>::iterator iter, iterl, iterp;
2644 iter = lista.begin();
2645 /*
2646 if(minipasma)
2647 {
2648 E=p0.en + (p1.en - p0.en)/256;
2649 }
2650 else
2651 {
2652 E=p0.en + (p1.en - p0.en)/2;
2653 }
2654 iter=lista.insert(++iter,punkt(E,czyosobliwa(E)));
2655 */
2656 double szer = dokl + 1; // zeby na wejsciu bylo wiecej niz dokl
2657 while(wynik.empty() && szer >= dokl)
2658 {
2659 iterp = lista.end();
2660 iterp--;
2661 while(iterp != lista.begin())
2662 {
2663 iterl = iterp;
2664 iterl--;
2665 E = (iterp->en + iterl->en) / 2;
2666 szer = iterp->en - iterl->en;
2667 // std::clog<<"El = "<<iterl->en<<" Eit = "<<E<<" Ep = "<<iterp->en<<"\n";
2668 iter = lista.insert(iterp, punkt(E, czyosobliwa(E)));
2669 if(znak * iter->wart < 0)
2670 {
2671 wynik.push_back(iterl->en);
2672 wynik.push_back(iter->en);
2673 wynik.push_back(iterp->en);
2674 break;
2675 }
2676 iterp = iterl;
2677 }
2678 }
2679 return wynik;
2680}
2681/*****************************************************************************/
2682 void struktura::profil(double Ek, double rozdz)
2683{
2684 double E0 = dol;
2685 if(Ek <= E0)
2686 {
2687 Error err; // MD
2688 err << "Zly zakres energii!\n";
2689 throw err;
2690 }
2691 for(double E = E0; E <= Ek; E += rozdz)
2692 {
2693 std::cout << E << "\t" << czyosobliwa(E) << "\n";
2694 }
2695 std::cout << std::flush;
2696}
2697/*****************************************************************************/
2698 void struktura::szukanie_poziomow(double Ek, double rozdz,
2699 bool debug) // Trzeba dorobic obsluge niezbiegania sie metody siecznych
2700{
2701 ofstream plikdebug, debugfun;
2702 if(debug)
2703 {
2704 plikdebug.open("debug");
2705 debugfun.open("debugfun");
2706 plikdebug << "E\tzera\n";
2707 plikdebug.precision(12);
2708 }
2709 double E0 = dol;
2710 if(Ek <= E0)
2711 {
2712 Error err; // MD
2713 err << "Zly zakres energii!\n";
2714 throw err;
2715 }
2716 int M = 2 * (kawalki.size() + 2) - 2;
2717 // Bylo 'int M = 2*kawalki.size() - 2;'
2718 double wartakt;
2719 double (struktura::*fun)(double) = &struktura::czyosobliwa;
2720 std::vector<double> trojka;
2721#ifdef LOGUJ // MD
2722 std::clog << "W szukaniu, E_pocz = " << Ek << "\n";
2723#endif
2724 double wartpop = czyosobliwa(Ek);
2725#ifdef LOGUJ // MD
2726 std::clog << "Pierwsza wartosc = " << wartpop << "\n";
2727#endif
2728 A2D V(M, M);
2729 int ilepoziomow;
2730 int liczbazer;
2731 int ostatnie_dobre, nast_dobre;
2732 double E;
2733 double ciagl, Eb0, Eb1;
2734 int licz; // licznik poprawek do E
2735 if(wartpop == 0)
2736 {
2737 Ek -= rozdz;
2738 wartpop = czyosobliwa(Ek);
2739 }
2740 E = Ek - rozdz;
2741 wartakt = czyosobliwa(E);
2742 while((E > E0) && (wartpop * wartakt > 0))
2743 {
2744 wartpop = wartakt;
2745 E -= rozdz;
2746 // std::clog<<"Ek = "<<E<<"\n";
2747 wartakt = czyosobliwa(E);
2748 }
2749 /* std::clog<<"Sieczne z krancami "<<E<<" i "<<(E + rozdz)<<"\n";
2750 double Eost = sieczne(fun, E, E + rozdz); // koniec szukania najwyzszego stanu */
2751#ifdef LOGUJ // MD
2752 std::clog << "Bisekcja z krancami " << E << " i " << (E + rozdz) << "\n";
2753#endif
2754 double Eost = bisekcja(fun, E, E + rozdz);
2755 if(debug)
2756 {
2757 plikdebug << Eost << "\t" << std::flush;
2758 }
2759 V = A2D(M, M);
2760 ilepoziomow = ilezer_ffal(Eost, V) + 1;
2761 if(debug)
2762 {
2763 plikdebug << (ilepoziomow - 1) << "\t" << (sprawdz_ciaglosc(Eost, V)) << std::endl;
2764 }
2765 ostatnie_dobre = ilepoziomow - 1;
2766 rozwiazania.resize(ilepoziomow);
2767 stan nowy(Eost, V, ostatnie_dobre);
2768 stan stymcz;
2769 rozwiazania[ilepoziomow - 1] = nowy;
2770#ifdef LOGUJ // MD
2771 std::clog << " Eost = " << Eost << " ilepoziomow = " << ilepoziomow << "\n";
2772#endif
2773 punkt ost, ostdob, pierw;
2774 nast_dobre = -1;
2775 while(ostatnie_dobre >= 1)
2776 {
2777 licz = 0;
2778 ostdob = punkt(rozwiazania[ostatnie_dobre]);
2779 E = ostdob.en - rozdz;
2780 ost = punkt(E, czyosobliwa(E));
2781 E = (nast_dobre >= 0) ? (rozwiazania[nast_dobre].poziom + rozdz) : (E0 + rozdz);
2782 pierw = punkt(E, czyosobliwa(E));
2783#ifdef LOGUJ // MD
2784 std::clog << "Eost = " << ost.en << " wartost = " << ost.wart << "\n";
2785 std::clog << "Epierw = " << pierw.en << " wartpierw = " << pierw.wart << "\n";
2786#endif
2787 if(ost.wart * pierw.wart > 0)
2788 {
2789#ifdef LOGUJ // MD
2790 std::clog << "Zageszczanie z pierw = " << pierw.en << " i ost = " << ost.en << "\n";
2791#endif
2792 trojka = zageszczanie(pierw, ost);
2793 for(int i = 1; i >= 0; i--) // bo sa dwa zera
2794 {
2795 // E = sieczne(fun, trojka[i], trojka[i+1]);
2796 E = bisekcja(fun, trojka[i], trojka[i + 1]);
2797 ciagl = sprawdz_ciaglosc(E, V);
2798 if(debug)
2799 {
2800 plikdebug << E << "\t" << ciagl << std::flush;
2801 }
2802 Eb0 = trojka[i];
2803 Eb1 = trojka[i + 1];
2804 while(ciagl < 0 || ciagl > 1e-7) // zmiana znaku na interfejsie albo marna ciaglosc
2805 {
2806 if(debug)
2807 {
2808 liczbazer = ilezer_ffal(E, V);
2809 stymcz = stan(E, V, liczbazer);
2810 funkcja1_do_pliku(debugfun, stymcz, 0.05);
2811 }
2812 licz++;
2813 if(licz > 10)
2814 {
2815 Error err; // MD
2816 err << "Dla E = " << E << " wychodzi niedobra funkcja, ciagl = " << ciagl << "\n";
2817 throw err;
2818 }
2819 if((this->*fun)(E) * (this->*fun)(Eb0) < 0)
2820 {
2821 Eb1 = E;
2822 }
2823 else
2824 {
2825 Eb0 = E;
2826 }
2827 E = bisekcja(fun, Eb0, Eb1);
2828 ciagl = sprawdz_ciaglosc(E, V);
2829 if(debug)
2830 {
2831 plikdebug << "poprawione E = " << E << " ciagl = " << ciagl << "\t" << std::flush;
2832 }
2833 }
2834 liczbazer = ilezer_ffal(E, V);
2835 if(debug)
2836 {
2837 plikdebug << "\t" << liczbazer << std::endl;
2838 }
2839#ifdef LOGUJ // MD
2840 std::clog << "E = " << E << "\tzer " << liczbazer << "\n";
2841#endif
2842 nowy = stan(E, V, liczbazer);
2843 if(liczbazer > ostatnie_dobre)
2844 {
2845 Error err; // MD
2846 err << "Za duzo zer!\n";
2847 throw err;
2848 }
2849 rozwiazania[liczbazer] = nowy;
2850 // if(debug && licz > 0) // rysuje funkcje jesli bylo podejrzane 0
2851 // funkcja1_do_pliku(debugfun, nowy, 0.1);
2852 }
2853 }
2854 else
2855 {
2856 // E = sieczne(fun, pierw.en, ost.en);
2857 E = bisekcja(fun, pierw.en, ost.en);
2858 ciagl = sprawdz_ciaglosc(E, V);
2859 if(debug)
2860 {
2861 plikdebug << E << "\t" << ciagl << std::flush;
2862 }
2863 Eb0 = pierw.en;
2864 Eb1 = ost.en;
2865 while(ciagl < 0 || ciagl > 1e-7) // zmiana znaku na interfejsie albo marna ciaglosc
2866 {
2867 if(debug)
2868 {
2869 liczbazer = ilezer_ffal(E, V);
2870 stymcz = stan(E, V, liczbazer);
2871 funkcja1_do_pliku(debugfun, stymcz, 0.05);
2872 }
2873 licz++;
2874 if(licz > 10)
2875 {
2876 Error err; // MD
2877 err << "Dla E = " << E << " wychodzi niedobra funkcja, ciagl = " << ciagl << "\n";
2878 throw err;
2879 }
2880 if((this->*fun)(E) * (this->*fun)(Eb0) < 0)
2881 {
2882 Eb1 = E;
2883 }
2884 else
2885 {
2886 Eb0 = E;
2887 }
2888 E = bisekcja(fun, Eb0, Eb1);
2889 ciagl = sprawdz_ciaglosc(E, V);
2890 if(debug)
2891 {
2892 plikdebug << "poprawione E = " << E << " ciagl = " << ciagl << "\t" << std::flush;
2893 }
2894 }
2895 liczbazer = ilezer_ffal(E, V);
2896 if(debug)
2897 {
2898 plikdebug << "\t" << liczbazer << std::endl;
2899 }
2900#ifdef LOGUJ // MD
2901 std::clog << "W else E = " << E << "\tzer " << liczbazer << "\n";
2902#endif
2903 nowy = stan(E, V, liczbazer);
2904 if(liczbazer > ostatnie_dobre)
2905 {
2906 Error err; // MD
2907 err << "Za duzo zer!\n";
2908 throw err;
2909 }
2910 if(rozwiazania[liczbazer].liczba_zer > -1) // juz jest stan z taka liczba zer
2911 {
2912 Error err; // MD
2913 err << "Bez sensu z zerami!\n";
2914 throw err;
2915
2916 /* stymcz = rozwiazania[liczbazer];
2917 if(stymcz.poziom > nowy.poziom)
2918 {
2919 if(nowy.liczba_zer == 0 || rozwiazania[liczbazer-1] > -1) // najnizszy poziom lub poniezej zajete
2920 {
2921 Error err; // MD
2922 err<<"Bardzo bez sensu z zerami!\n";
2923 throw err;
2924 }
2925 nowy.liczba_zer -= 1;
2926 rozwiazania[nowy.liczba_zer] = nowy;
2927 // std::clog<<"Poziom o energii
2928 }
2929 else
2930 {
2931 if(stymcz.liczba_zer == 0 || rozwiazania[liczbazer-1] > -1)
2932 {
2933 Error err; // MD
2934 err<<"Bardzo bez sensu z zerami!\n";
2935 throw err;
2936 }
2937 stymcz.liczba_zer -= 1;
2938 rozwiazania[stymcz.liczba_zer] = stymcz;
2939 rozwiazania[nowy.liczba_zer] = nowy;
2940 }
2941 */
2942 }
2943 else
2944 {
2945 rozwiazania[liczbazer] = nowy;
2946 }
2947 if(debug && licz > 0) // rysuje funkcje jesli bylo podejrzane 0
2948 funkcja1_do_pliku(debugfun, nowy, 0.1);
2949 }
2950#ifdef LOGUJ // MD
2951 std::clog << "ostatnie_dobre = " << ostatnie_dobre << "\n";
2952#endif
2953 while(ostatnie_dobre >= 1 && rozwiazania[ostatnie_dobre - 1].liczba_zer >= 0)
2954 {
2955 ostatnie_dobre--;
2956 }
2957 nast_dobre = ostatnie_dobre - 1;
2958 while(nast_dobre >= 0 && rozwiazania[nast_dobre].liczba_zer < 0)
2959 {
2960 nast_dobre--;
2961 }
2962#ifdef LOGUJ // MD
2963 std::clog << "ostatnie_dobre = " << ostatnie_dobre << " nastepne_dobre = " << nast_dobre << "\n";
2964#endif
2965 }
2966 plikdebug.close();
2967 debugfun.close();
2968#ifdef LOGUJ // MD
2969 std::clog << "Liczba rozwiazan = " << rozwiazania.size() << "\n";
2970#endif
2971}
2972/*****************************************************************************/
2973 void struktura::szukanie_poziomow_2(double Ek, double rozdz, bool debug)
2974{
2975 ofstream plikdebug, debugfun;
2976 if(debug)
2977 {
2978 plikdebug.open("debug");
2979 debugfun.open("debugfun");
2980 plikdebug << "E\tzera\n";
2981 plikdebug.precision(12);
2982 }
2983 std::list<stan> listast; // lista znalezionych stanow
2984 std::list<stan>::iterator itostd,
2985 itnastd; // najwyzszy stan nad kto(c)ym wszystko jest znalezione, nastepny stan, nad ktorym trzeba jeszcze szukac
2986 std::list<stan>::iterator iter; // pomocniczy
2987 double E0 = dol;
2988 if(Ek <= E0)
2989 {
2990 Error err; // MD
2991 err << "Zly zakres energii!\n";
2992 throw err;
2993 }
2994 int M = 2 * (kawalki.size() + 2) - 2;
2995 // Bylo 'int M = 2*kawalki.size() - 2;'
2996 double wartakt;
2997 double (struktura::*fun)(double) = &struktura::czyosobliwa;
2998 std::vector<double> trojka;
2999#ifdef LOGUJ // MD
3000 std::clog << "W szukaniu, E_pocz = " << Ek << "\n";
3001#endif
3002 double wartpop = czyosobliwa(Ek);
3003#ifdef LOGUJ // MD
3004 std::clog << "Pierwsza wartosc = " << wartpop << "\n";
3005#endif
3006 A2D V(M, M);
3007 //1.5
3008 std::vector<parad> wspolcz;
3009 std::ofstream pliktestowy;//("funkcjfalowe");
3010 //
3011 int ilepoziomow;
3012 int liczbazer;
3013 // int ostatnie_dobre , nast_dobre;
3014 double E;
3015 double B0, B1;
3016 // double ciagl, Eb0, Eb1;
3017 // int licz; // licznik poprawek do E
3018 if(wartpop == 0)
3019 {
3020 Ek -= rozdz;
3021 wartpop = czyosobliwa(Ek);
3022 }
3023 E = Ek - rozdz;
3024 wartakt = czyosobliwa(E);
3025 while((E > E0) && (wartpop * wartakt > 0))
3026 {
3027 wartpop = wartakt;
3028 //1.5
3029 // E -= rozdz;
3030 E -= rozdz*1e3; // żeby nie czekać wieki
3031 //
3032 // E -= 10*rozdz;
3033 // std::clog<<"Ek = "<<E<<"\n";
3034 wartakt = czyosobliwa(E);
3035 }
3036 /* std::clog<<"Sieczne z krancami "<<E<<" i "<<(E + rozdz)<<"\n";
3037 double Eost = sieczne(fun, E, E + rozdz); // koniec szukania najwyzszego stanu */
3038 std::clog<<"Bisekcja z krancami "<<E<<" i "<<(E + rozdz*1e3)<<"\n";
3039 double Eost = bisekcja(fun, E, E + rozdz*1e3);
3040
3041 // dopiski 1.5:
3042 parad ABkon = sklejanie_od_lewej(Eost);
3043 // std::clog<<"NOWE: A, B = "<<ABkon.first<<" "<<ABkon.second<<std::endl;
3044 std::clog<<"STARE: A, B = "<<ABkon.first<<" "<<ABkon.second<<"\tE = "<<Eost<<std::endl;
3045 B0 = ABkon.second;
3046 double nowaE = poprawianie_poziomu(Eost, 1e-11);
3047 ABkon = sklejanie_od_lewej(nowaE);
3048 B1 = ABkon.second;
3049 std::clog<<"NOWSZE: A, B = "<<ABkon.first<<" "<<ABkon.second<<"\tE = "<<nowaE<<std::endl;
3050 if(abs(B1) > abs(B0))
3051 {
3052 std::clog<<"B sie pogorszyło(0)! B0, B1 "<<B0<<", "<<B1<<std::endl;
3053 }
3054 std::clog<<"zmiana E = "<<(Eost-nowaE)<<std::endl;
3055 Eost = nowaE;
3056 //
3057
3058 if(debug)
3059 {
3060 plikdebug << Eost << "\t" << std::flush;
3061 }
3062 V = A2D(M, M);
3063 ilepoziomow = ilezer_ffal(Eost, V) + 1;
3064 wspolcz = wsp_sklejanie_od_lewej(Eost);
3065 int nowezera = ilezer_ffal(Eost, wspolcz);
3066 std::clog<<"stara i nowa liczba zer "<<ilepoziomow-1<<", "<<nowezera<<std::endl;
3067 if(debug)
3068 {
3069 plikdebug << (ilepoziomow - 1) << "\t" << (sprawdz_ciaglosc(Eost, V)) << std::endl;
3070 }
3071 // ostatnie_dobre = ilepoziomow - 1;
3072 stan nowy(Eost, V, ilepoziomow - 1);
3073 //1.5
3074 stan nowszy = stan(Eost, wspolcz, ilepoziomow - 1);
3075 std::string temat = "funkcjafalowa_E";
3076 std::string nazwapliku = temat + std::to_string(ilepoziomow-1);
3077 pliktestowy.open(nazwapliku);
3078 funkcja1_do_pliku(pliktestowy, nowy, 1e-2); //
3079 pliktestowy.close();
3080 pliktestowy.open(nazwapliku + "prim");
3081 funkcja1_do_pliku(pliktestowy, nowszy, 1e-2); //
3082 pliktestowy.close();
3083 pliktestowy.open("studnie");
3084 struktura_do_pliku(pliktestowy);
3085 pliktestowy.close();
3086 //
3087 listast.push_back(nowy);
3088 itostd = listast.end();
3089 --itostd; // iterator na ostatni element
3090 stan stymcz;
3091#ifdef LOGUJ // MD
3092 std::clog << " Eost = " << Eost << " ilepoziomow = " << ilepoziomow << "\n";
3093#endif
3094 punkt ost, ostdob, pierw;
3095 itnastd = listast.begin();
3096 // while(listast.size() < ilepoziomow)
3097 // while(itostd != listast.begin())
3098 while((int)listast.size() < ilepoziomow)
3099 {
3100 // licz = 0;
3101 ostdob = *itostd;
3102 E = ostdob.en - rozdz;
3103 ost = punkt(E, czyosobliwa(E));
3104 // E = (nast_dobre >= 0)?(rozwiazania[nast_dobre].poziom + rozdz):(E0 + rozdz);
3105 E = (itnastd != itostd) ? (itnastd->poziom + rozdz) : (E0 + rozdz);
3106 pierw = punkt(E, czyosobliwa(E));
3107
3108#ifdef LOGUJ // MD
3109 std::clog << "Eost = " << ost.en << " wartost = " << ost.wart << "\n";
3110 std::clog << "Epierw = " << pierw.en << " wartpierw = " << pierw.wart << "\n";
3111#endif
3112
3113 if(ost.wart * pierw.wart > 0)
3114 {
3115#ifdef LOGUJ // MD
3116 std::clog << "Zageszczanie z pierw = " << pierw.en << " i ost = " << ost.en << "\n";
3117#endif
3118 trojka = zageszczanie(pierw, ost);
3119 if(!trojka.empty())
3120 {
3121 iter = itostd;
3122 for(int i = 1; i >= 0; i--) // bo sa dwa zera
3123 {
3124 // E = sieczne(fun, trojka[i], trojka[i+1]);
3125 E = bisekcja(fun, trojka[i], trojka[i + 1]);
3126
3127 //1.5
3128 /* ABkon = sklejanie_od_lewej(E);
3129 B0 = ABkon.second;
3130 std::clog<<"STARE: A, B = "<<ABkon.first<<" "<<ABkon.second<<std::endl;*/
3131 nowaE = poprawianie_poziomu(E, 1e-11);
3132 /* ABkon = sklejanie_od_lewej(nowaE);
3133 B1 = ABkon.second;
3134 std::clog<<"NOWSZE: A, B = "<<ABkon.first<<" "<<ABkon.second<<"\tE = "<<nowaE<<std::endl;
3135 if(abs(B1) > abs(B0))
3136 {
3137 std::clog<<"B sie pogorszyło(1)! B0, B1 "<<B0<<", "<<B1<<std::endl;
3138 }
3139
3140 std::clog<<"zmiana E = "<<(E-nowaE)<<std::endl; */
3141
3142 E = nowaE;
3143 //
3144
3145 liczbazer = ilezer_ffal(E, V);
3146#ifdef LOGUJ // MD
3147 std::clog << "E = " << E << " zer jest " << liczbazer << "\n";
3148#endif
3149 nowy = stan(E, V, liczbazer);
3150
3151 //1.5
3152 wspolcz = wsp_sklejanie_od_lewej(E);
3153 nowezera = ilezer_ffal(E, wspolcz);
3154 std::clog<<"\nstara i nowa liczba zer "<<liczbazer<<", "<<nowezera<<std::endl;
3155
3156 nazwapliku = temat + std::to_string(liczbazer);
3157 pliktestowy.open(nazwapliku);
3158 funkcja1_do_pliku(pliktestowy, nowy, 1e-2);
3159 pliktestowy.close();
3160 nowszy = stan(E, wspolcz, nowezera);
3161 pliktestowy.open(nazwapliku + "prim");
3162 funkcja1_do_pliku(pliktestowy, nowszy, 1e-2); //
3163 pliktestowy.close();
3164
3165 //
3166 // funkcja1_do_pliku(pliktestowy, nowy, 1e-2);
3167 /* if(liczbazer > ostatnie_dobre)
3168 {
3169 Error err; // MD
3170 err<<"Za duzo zer!\n";
3171 throw err;
3172 }
3173 */
3174 listast.insert(iter, nowy);
3175 --iter;
3176 // rozwiazania[liczbazer] = nowy;
3177 // if(debug && licz > 0) // rysuje funkcje jesli bylo podejrzane 0
3178 // funkcja1_do_pliku(debugfun, nowy, 0.1);
3179 }
3180 }
3181 else
3182 {
3183 itostd = itnastd;
3184#ifdef LOGUJ // MD
3185 std::clog << "Nie znalazl a powinien\n";
3186#endif
3187 }
3188 }
3189 else
3190 {
3191 // E = sieczne(fun, pierw.en, ost.en);
3192 E = bisekcja(fun, pierw.en, ost.en);
3193
3194 //1.5
3195 ABkon = sklejanie_od_lewej(E);
3196
3197 B0 = ABkon.second;
3198 std::clog<<"STARE: A, B = "<<ABkon.first<<" "<<ABkon.second<<"\tE = "<<E<<std::endl;
3199 nowaE = poprawianie_poziomu(E, 1e-11);
3200 ABkon = sklejanie_od_lewej(nowaE);
3201 B1 = ABkon.second;
3202 std::clog<<"NOWSZE: A, B = "<<ABkon.first<<" "<<ABkon.second<<"\tE = "<<nowaE<<std::endl;
3203
3204 if(abs(B1) > abs(B0))
3205 {
3206 std::clog<<"B sie pogorszyło(2)! B0, B1 "<<B0<<", "<<B1<<std::endl;
3207 }
3208
3209
3210 std::clog<<"zmiana E = "<<(E-nowaE)<<std::endl;
3211 E = nowaE;
3212 //
3213
3214 liczbazer = ilezer_ffal(E, V);
3215#ifdef LOGUJ // MD
3216 std::clog << "E = " << E << " zer jest " << liczbazer << "\n";
3217#endif
3218 nowy = stan(E, V, liczbazer);
3219 // funkcja1_do_pliku(pliktestowy, nowy, 1e-2);
3220 if(liczbazer >= itostd->liczba_zer)
3221 {
3222 Error err; // MD
3223 err << "Za duzo zer!\n";
3224 throw err;
3225 }
3226 listast.insert(itostd, nowy);
3227 // rozwiazania[liczbazer] = nowy;
3228 }
3229 // std::clog<<"ostatnie_dobre = "<<ostatnie_dobre<<"\n";
3230 // while(ostatnie_dobre >= 1 && rozwiazania[ostatnie_dobre - 1].liczba_zer >= 0 )
3231 while((itostd != listast.begin()) &&
3232 (itostd->liczba_zer <= std::prev(itostd)->liczba_zer + 1)) // stany na liscie schodza po jeden w liczbie zer
3233 {
3234#ifdef LOGUJ // MD
3235 if(itostd->liczba_zer < std::prev(itostd)->liczba_zer + 1)
3236 {
3237 std::clog << "\nKlopot z miejscami zerowymi: E1 = " << (itostd->poziom) << " zer " << (itostd->liczba_zer)
3238 << "\nE2 = " << (std::prev(itostd)->poziom) << " zer " << (std::prev(itostd)->liczba_zer)
3239 << "\n";
3240 }
3241#endif
3242 --itostd;
3243 }
3244 // kontrolne wypisywanie
3245 /* iter = listast.end();
3246 while(iter != listast.begin())
3247 {
3248 --iter;
3249 // std::clog<<"zera = "<<(iter->liczba_zer)<<"\tE = "<<(iter->poziom)<<"\n";
3250 }
3251 // koniec kontrolnego wyppisywania
3252 */
3253 itnastd = (itostd == listast.begin()) ? listast.begin() : std::prev(itostd);
3254#ifdef LOGUJ // MD
3255 std::clog << "ostatnie_dobre = " << (itostd->poziom) << " nastepne_dobre = " << (itnastd->poziom) << "\n";
3256#endif
3257 }
3258 if(debug)
3259 {
3260 plikdebug.close();
3261 debugfun.close();
3262 }
3263#ifdef LOGUJ // MD
3264 std::clog << "Liczba rozwiazan = " << rozwiazania.size() << "\n";
3265#endif
3266 rozwiazania.reserve(listast.size());
3267 iter = listast.begin();
3268 while(iter != listast.end())
3269 {
3270 rozwiazania.push_back(*iter);
3271 // std::clog<<"zera = "<<(iter->liczba_zer)<<"\tE = "<<(iter->poziom)<<"\n";
3272 ++iter;
3273 }
3274}
3275/*****************************************************************************/
3276 void struktura::stany_z_tablicy(const std::vector<double> & energie)
3277{
3278 int M = 2 * (kawalki.size() + 2) - 2;
3279 A2D V(M, M);
3280 double liczbazer;
3281 double E;
3282 stan nowy;
3283 rozwiazania.reserve(energie.size());
3284 for(int i = 0; i <= (int)energie.size() - 1; ++i)
3285 {
3286 E = energie[i];
3287 liczbazer = ilezer_ffal(E, V);
3288#ifdef LOGUJ // MD
3289 if(liczbazer != i)
3290 {
3291 std::cerr << "E = " << E << " zer jest " << liczbazer << " a powinno byc " << i << "\n";
3292 // throw err;
3293 }
3294#endif
3295 nowy = stan(E, V, liczbazer);
3296 rozwiazania.push_back(nowy);
3297 }
3298}
3299/*****************************************************************************/
3300 double struktura::sieczne(double (struktura::*f)(double), double pocz, double kon)
3301{
3302#ifdef LOGUJ // MD
3303 std::clog.precision(12);
3304#endif
3305 // const double eps = 1e-14; // limit zmian x
3306 double dokl = 1e-10;
3307 double xp = kon;
3308 double xl = pocz;
3309#ifdef LOGUJ // MD
3310 std::clog << "Sieczne. Wartosci na koncach: " << ((this->*f)(pocz)) << ", " << ((this->*f)(kon)) << "\n";
3311#endif
3312 if((this->*f)(pocz) * (this->*f)(kon) > 0)
3313 {
3314 Error err; // MD
3315 err << "Zle krance przedzialu!\n";
3316 throw err;
3317 }
3318
3319 double x, fc, fp, fl, xlp, xpp; // -p -- poprzednie krance
3320 double fpw, flw; // krance funkcji do wzoru -- moga roznic sie od prawdziwych w Illinois
3321 fl = (this->*f)(xl);
3322 fp = (this->*f)(xp);
3323 flw = fl;
3324 fpw = fp;
3325 xlp = (xl + xp) / 2;
3326 xpp = xlp; // zeby na pewno bylo na poczatku rozne od prawdzwych koncow
3327 do
3328 {
3329 x = xp - fpw * (xp - xl) / (fpw - flw);
3330 fpw = fp;
3331 flw = fl;
3332 fc = (this->*f)(x);
3333 if(fc == 0)
3334 {
3335 break;
3336 }
3337 if(fc * fl < 0) // c bedzie nowym prawym koncem
3338 {
3339 // std::clog<<"xlp - xl = "<<(xlp - xl)<<"\n";
3340 if(xlp == xl) // trzeba poprawic ten kraniec (metoda Illinois)
3341 {
3342#ifdef LOGUJ // MD
3343 std::clog << "Lewy Illinois\n";
3344#endif
3345 flw = flw / 2;
3346 }
3347 xpp = xp;
3348 xp = x;
3349 xlp = xl;
3350 fp = fc;
3351 }
3352 else // c bedzie nowym lewym koncem
3353 {
3354 // std::clog<<"xpp - xp = "<<(xpp - xp)<<"\n";
3355 if(xpp == xp) // trzeba poprawic ten kraniec (metoda Illinois)
3356 {
3357#ifdef LOGUJ // MD
3358 std::clog << "Prawy Illinois\n";
3359#endif
3360 fpw = fpw / 2;
3361 }
3362 xlp = xl;
3363 xl = x;
3364 xpp = xp;
3365 fl = fc;
3366 // fp = (this->*f)(kon);
3367 }
3368#ifdef LOGUJ // MD
3369 std::clog << "x = " << x << "\tf(x) = " << fc << "\txl = " << xl << " xp = " << xp << " f(xl) = " << fl
3370 << " f(xp) = " << fp << "\n"; //<<"\txp - x = "<<(xp - x)<<"\n";
3371#endif
3372 }
3373 while(xp - xl >= dokl);
3374#ifdef LOGUJ // MD
3375 if(fc * fc > 1e-8)
3376 std::cerr << "\nfc = " << fc << " zamiast 0!\n";
3377#endif
3378 return x;
3379}
3380/*****************************************************************************/
3381//1.5
3382//double struktura::bisekcja(double (struktura::*f)(double), double pocz, double kon,)
3383 double struktura::bisekcja(double (struktura::*f)(double), double pocz, double kon, double dokl) // domyślnie dokl = 1e-9
3384//
3385{
3386#ifdef LOGUJ // MD
3387 std::clog.precision(12);
3388#endif
3389 // const double eps = 1e-14; // limit zmian x
3390 //1.5
3391 // double dokl = 1e-9;
3392 //
3393 // int limpetli = 30;
3394 double xp = kon;
3395 double xl = pocz;
3396 // std::clog<<"Bisekcja. Wartości na końcach: "<<((this->*f)(pocz))<<", "<<((this->*f)(kon))<<"\tdokl = "<<dokl<<"\n";
3397 if((this->*f)(pocz) * (this->*f)(kon) > 0)
3398 {
3399 Error err; // MD
3400 err << "Zle krance przedzialu!\n";
3401 throw err;
3402 }
3403
3404 if(pocz >= kon)
3405 {
3406 std::cerr<<"Zła kolejność krańców!\n";
3407 abort();
3408 }
3409 //
3410
3411 double x, fc, fl; // -p -- poprzednie krance
3412 //1.5
3413 double fp;
3414 //
3415 fl = (this->*f)(xl);
3416 int liczpetli = 0;
3417 // fp = (this->*f)(xp); //Wystarczy chyba tylko tu, ale wtedy nie wypisze dobrze wartosci na krancach
3418 fp = (this->*f)(xp); //Wystarczy chyba tylko tu, ale wtedy nie wypisze dobrze wartości na krańcach
3419 //
3420 do
3421 {
3422 x = (xp + xl) / 2;
3423 fc = (this->*f)(x);
3424 if(fc == 0)
3425 {
3426 break;
3427 }
3428 if(fc * fl < 0) // c bedzie nowym prawym koncem
3429 {
3430 xp = x;
3431 // fp = (this->*f)(xp);
3432 //1.5
3433 fp = fc;
3434 //
3435 }
3436 else // c bedzie nowym lewym koncem
3437 {
3438 xl = x;
3439 //1.5
3440 // fl = (this->*f)(xl);
3441 fl = fc;
3442 //
3443 }
3444 liczpetli += 1;
3445 // std::clog<<"petla "<<liczpetli<<" x = "<<x<<"\tf(x) = "<<fc<<"\txl = "<<xl<<" xp = "<<xp<<"\n"; //f(xl) =
3446 // "<<fl<<" f(xp) = "<<fp<<"\n";//<<"\txp - x = "<<(xp - x)<<"\n";
3447 }
3448 // while(xp - xl >= dokl);
3449 // while(liczpetli < limpetli && (xp - xl >= dokl || fc*fc > 1e-8) );
3450 // while(liczpetli < limpetli && xp - xl >= dokl );
3451 while(xp - xl >= dokl);
3452 /*
3453 if(liczpetli >= limpetli)
3454 {
3455 std::cerr<<"\nPrzekroczony limit liczby petli\nOsiagnieta dokladnosc:"<<(xp - xl)<<"\n";
3456 }
3457 */
3458 //1.5
3459 // if(fc*fc > 1e-8)
3460 // {
3461 // std::cerr<<"\nfc = "<<fc<<" zamiast 0!\n";
3462 // }
3463
3464 // return x;
3465 // std::cerr<<"\nfl = "<<fl<<" fp = "<<fp<<"\n";
3466 // std::cerr<<"Dx = "<<xp-xl<<"\n";
3467 // return (abs(fl) > abs(fp))?xp:xl;
3468 //usiecznienie
3469 double wsp = (fp -fl)/(xp - xl);
3470 double xsiecz = xl - fl/wsp;
3471 double wyn;
3472 double fsiecz = (this->*f)(xsiecz);
3473 if(abs(fsiecz) < abs(fp) && abs(fsiecz) < abs(fl)) // z siecznej lepsze
3474 wyn = xsiecz;
3475 else
3476 {
3477 if(abs(fp) > abs(fl)) // który krzniec lepszy?
3478 wyn = xl;
3479 else
3480 wyn = xp;
3481 }
3482 return wyn;
3483 //
3484}
3485/*****************************************************************************/
3486 double struktura::norma_stanu(stan & st) // liczy norme i wypelnia prawdopodobienstwa
3487{
3488 double porcja = lewa.norma_kwadr(st.poziom, st.wspolczynniki.front());
3489 st.prawdopodobienstwa.push_back(porcja);
3490 double norma2 = porcja;
3491 int i; // LUKASZ Visual 6
3492 for(/*int*/ i = 0; i <= (int)kawalki.size() - 1; i++) // LUKASZ Visual 6
3493 {
3494 porcja = kawalki[i].norma_kwadr(st.poziom, st.wspolczynniki[2 * i + 1], st.wspolczynniki[2 * i + 2]);
3495 st.prawdopodobienstwa.push_back(porcja);
3496 norma2 += porcja;
3497 }
3498 porcja = prawa.norma_kwadr(st.poziom, st.wspolczynniki.back());
3499 st.prawdopodobienstwa.push_back(porcja);
3500 norma2 += porcja;
3501 for(/*int*/ i = 0; i <= (int)st.prawdopodobienstwa.size() - 1; i++) // LUKASZ Visual 6
3502 {
3503 st.prawdopodobienstwa[i] /= norma2;
3504 }
3505 return sqrt(norma2);
3506}
3507/*****************************************************************************/
3508 void struktura::normowanie()
3509{
3510 std::vector<stan>::iterator it = rozwiazania.begin();
3511 // std::clog<<"Liczba stanow = "<<rozwiazania.size()<<"\n";
3512 double norma;
3513 while(it != rozwiazania.end())
3514 {
3515 norma = norma_stanu(*it);
3516 // std::clog<<"Norma dla E = "<<(it->poziom)<<" wynosi "<<norma<<"\n";
3517 for(int i = 0; i <= (int)it->wspolczynniki.size() - 1; i++)
3518 {
3519 it->wspolczynniki[i] /= norma;
3520 }
3521 it++;
3522 }
3523}
3524/*****************************************************************************/
3525 double struktura::ilenosnikow(double qFl, double T)
3526{
3527 // std::cout << "\tjestem w ilenosnikow (2 argumenty): " << "\n"; // TEST // TEST LUKASZ
3528 double tylenosnikow = 0.0;
3529 double niepomnozone = 0.0;
3530 double calkazFD;
3531 double szer;
3532 double gam32 = sqrt(pi) / 2; // Gamma(3/2)
3533 double kT = kB * T;
3534 std::vector<stan>::iterator it = rozwiazania.end();
3535 // std::clog<<"Liczba stanow = "<<rozwiazania.size()<<"\n";
3536 // std::clog<<"rodzaj = "<<typ<<"\n";
3537 while(it != rozwiazania.begin()) // suma po poziomach
3538 {
3539 --it;
3540 // std::clog<<"niepomonzone = "<<niepomnozone<<"\n";
3541 calkazFD = kT * log(1 + exp((qFl - it->poziom) / kT));
3542 niepomnozone = lewa.norma_kwadr(it->poziom, it->wspolczynniki.front()) * lewa.masa_r;
3543 niepomnozone += prawa.norma_kwadr(it->poziom, it->wspolczynniki.back()) * prawa.masa_r;
3544 for(int i = 0; i <= (int)kawalki.size() - 1; i++) // Suma po warstwach
3545 {
3546 // std::clog<<"niepomonzone = "<<niepomnozone<<"\n";
3547 niepomnozone +=
3548 kawalki[i].norma_kwadr(it->poziom, it->wspolczynniki[2 * i + 1], it->wspolczynniki[2 * i + 2]) *
3549 kawalki[i].masa_r;
3550 }
3551 tylenosnikow += niepomnozone * calkazFD / pi; // spiny juz sa
3552 }
3553 niepomnozone = 0;
3554 double Egorna = lewa.y;
3555 for(int i = 0; i <= (int)kawalki.size() - 1; i++) // Po stanach 3D
3556 {
3557 szer = kawalki[i].x_kon - kawalki[i].x_pocz;
3558 niepomnozone += szer * sqrt(2 * kawalki[i].masa_p(Egorna)) * kawalki[i].masa_r; // Spin juz jest
3559 }
3560 tylenosnikow += niepomnozone * kT * gam32 * sqrt(kT) * 2 * gsl_sf_fermi_dirac_half((qFl - Egorna) / (kB * T)) /
3561 (2 * pi * pi); // Powinno byc rozpoznawanie, ktore warstwy wystaja ponad poziom bariery, ale w GSL nie ma
3562 // niezupelnych calek F-D. Mozna by je przyblizyc niezupelnymi funkcjami Gamma.
3563 return tylenosnikow;
3564}
3565/*****************************************************************************/
3566 double struktura::ilenosnikow(double qFl, double T, std::set<int> ktore_warstwy)
3567{
3568 // std::cout << "\tjestem w ilenosnikow (3 argumenty): " << "\n"; // TEST // TEST LUKASZ
3569 double tylenosnikow = 0.0;
3570 double niepomnozone = 0.0;
3571 double calkazFD;
3572 double szer;
3573 double gam32 = sqrt(pi) / 2; // Gamma(3/2)
3574 double kT = kB * T;
3575 std::vector<stan>::iterator it = rozwiazania.end();
3576 // std::clog<<"Liczba stanow = "<<rozwiazania.size()<<"\n";
3577 std::set<int>::iterator setit;
3578 while(it != rozwiazania.begin()) // suma po poziomach
3579 {
3580 --it;
3581 calkazFD = kT * log(1 + exp((qFl - it->poziom) / kT));
3582 niepomnozone = 0.0;
3583 /* niepomnozone = lewa.norma_kwadr(it->poziom, it->wspolczynniki.front()) * lewa.masa_r;
3584 niepomnozone += prawa.norma_kwadr(it->poziom, it->wspolczynniki.back()) * prawa.masa_r;
3585 */
3586 for(setit = ktore_warstwy.begin(); setit != ktore_warstwy.end(); ++setit) // Suma po warstwach
3587 {
3588 // std::clog<<"nosniki w warstwie "<<(*setit)<<"\n";
3589 // std::cout << "\tit->poziom: " << it->poziom << "\n"; // TEST // TEST LUKASZ
3590 // std::cout << "\tit->wspolczynniki[2*(*setit) + 1]: " << it->wspolczynniki[2*(*setit) + 1] << "\n"; // TEST
3591 // // TEST LUKASZ std::cout << "\tit->wspolczynniki[2*(*setit) + 2]: " << it->wspolczynniki[2*(*setit) + 2] <<
3592 // "\n"; // TEST // TEST LUKASZ if (abs(it->wspolczynniki[2*(*setit) + 2]) < 1e-6) continue; // DODALEM LINIE
3593 // LUKASZ std::cout << "\tkawalki[*setit].masa_r: " << kawalki[*setit].masa_r << "\n"; // TEST // TEST LUKASZ
3594 // if (!kawalki[*setit].masa_r) continue; // DODALEM LINIE LUKASZ
3595 // std::cout << "\tkawalki[*setit].masa_r_po: " << kawalki[*setit].masa_r << "\n"; // TEST // TEST LUKASZ
3596 niepomnozone += kawalki[*setit].norma_kwadr(it->poziom, it->wspolczynniki[2 * (*setit) + 1],
3597 it->wspolczynniki[2 * (*setit) + 2]) *
3598 kawalki[*setit].masa_r;
3599 // std::cout << "\tniepomnozone: " << niepomnozone << "\n"; // TEST // TEST LUKASZ
3600 }
3601 tylenosnikow += niepomnozone * calkazFD / pi; // spiny juz sa
3602 }
3603 niepomnozone = 0;
3604 double Egorna = lewa.y;
3605 for(setit = ktore_warstwy.begin(); setit != ktore_warstwy.end(); ++setit)
3606 {
3607 szer = kawalki[*setit].x_kon - kawalki[*setit].x_pocz;
3608 niepomnozone += szer * sqrt(2 * kawalki[*setit].masa_p(Egorna)) * kawalki[*setit].masa_r; // Spin juz jest
3609 }
3610 tylenosnikow += niepomnozone * kT * gam32 * sqrt(kT) * 2 * gsl_sf_fermi_dirac_half((qFl - Egorna) / kT) /
3611 (2 * pi * pi); // Powinno byc rozpoznawanie, ktore warstwy wystaja ponad poziom bariery, ale w GSL nie ma
3612 // niezupelnych calek F-D. Mozna by je przyblizyc niezupelnymi funkcjami Gamma.
3613 return tylenosnikow;
3614}
3615/*****************************************************************************/
3616std::vector<double> struktura::koncentracje_w_warstwach(double qFl, double T)
3617{
3618 double tylenosnikow = 0;
3619 double niepomnozone;
3620 double calkazFD, calkaFD12;
3621 double szer;
3622 double gam32 = sqrt(pi) / 2; // Gamma(3/2)
3623 double kT = kB * T;
3624 double Egorna = lewa.y;
3625 calkaFD12 = gsl_sf_fermi_dirac_half((qFl - Egorna) / (kB * T));
3626 std::vector<double> koncentr(kawalki.size() + 2);
3627 std::vector<stan>::iterator it;
3628 // std::clog<<"Liczba stanow = "<<rozwiazania.size()<<"\n";
3629 koncentr[0] = sqrt(2 * lewa.masa_p) * lewa.masa_r * kT * gam32 * sqrt(kT) * 2 * calkaFD12 / (2 * pi * pi);
3630 for(int i = 0; i <= (int)kawalki.size() - 1; i++) // Suma po warstwach
3631 {
3632 it = rozwiazania.end();
3633 // std::clog<<"i = "<<i<<"\n";
3634 tylenosnikow = 0;
3635 while(it != rozwiazania.begin()) // suma po poziomach
3636 {
3637 --it;
3638 if(i == 1)
3639 {
3640 // std::clog<<"Zawartosc dla poziomu "<<(it->poziom)<<" wynosi
3641 //"<<(kawalki[i].norma_kwadr(it->poziom, it->wspolczynniki[2*i + 1], it->wspolczynniki[2*i + 2]))<<"\n";
3642 }
3643 calkazFD = kT * log(1 + exp((qFl - it->poziom) / kT));
3644 niepomnozone =
3645 kawalki[i].norma_kwadr(it->poziom, it->wspolczynniki[2 * i + 1], it->wspolczynniki[2 * i + 2]) *
3646 kawalki[i].masa_r;
3647 tylenosnikow += niepomnozone * calkazFD / pi; // spiny juz sa
3648 }
3649 szer = kawalki[i].x_kon - kawalki[i].x_pocz;
3650 koncentr[i + 1] = tylenosnikow / szer +
3651 sqrt(2 * kawalki[i].masa_p(Egorna)) * kawalki[i].masa_r * kT * gam32 * sqrt(kT) * 2 * calkaFD12 / (2 * pi * pi);
3652 }
3653 koncentr.back() = koncentr.front();
3654 return koncentr;
3655}
3656/*****************************************************************************/
3657 void struktura::struktura_do_pliku(std::ofstream & plik)
3658{
3659 double X = 0.1; // do zamiany A->nm // dodalem X, LUKASZ
3660 std::vector<warstwa>::iterator it_war = kawalki.begin();
3661 plik << dlugosc_na_A(lewa.iks) << " " << (lewa.y) << "\n";
3662 while(it_war != kawalki.end())
3663 {
3664 plik << dlugosc_na_A(it_war->x_pocz) << " " << (it_war->y_pocz) << "\n";
3665 plik << dlugosc_na_A(it_war->x_kon) << " " << (it_war->y_kon) << "\n";
3666 it_war++;
3667 }
3668 plik << dlugosc_na_A(prawa.iks) << " " << (prawa.y);
3669}
3670/*****************************************************************************/
3671 void struktura::poziomy_do_pliku_(std::ofstream & plik, char c, double iRefBand, double iX1,
3672 double iX2) // LUKASZ dodalem cala funkcje
3673{
3674 std::vector<stan>::iterator it_stan = rozwiazania.begin();
3675 plik << iX1 * 0.1;
3676 while(it_stan != rozwiazania.end())
3677 {
3678 if(c == 'e')
3679 plik << '\t' << iRefBand + (it_stan->poziom); // TEST1
3680 else if(c == 'h')
3681 plik << '\t' << iRefBand - (it_stan->poziom); // TEST1
3682 // if (c == 'e') plik << '\t' << 0*iRefBand + (it_stan->poziom); //TEST2
3683 // else if (c == 'h') plik << '\t' << 0*iRefBand - (it_stan->poziom); //TEST2
3684 it_stan++;
3685 }
3686 plik << '\n';
3687
3688 it_stan = rozwiazania.begin();
3689 plik << iX2 * 0.1;
3690 while(it_stan != rozwiazania.end())
3691 {
3692 if(c == 'e')
3693 plik << '\t' << iRefBand + (it_stan->poziom); // TEST1
3694 else if(c == 'h')
3695 plik << '\t' << iRefBand - (it_stan->poziom); // TEST1
3696 // if (c == 'e') plik << '\t' << 0*iRefBand + (it_stan->poziom); //TEST2
3697 // else if (c == 'h') plik << '\t' << 0*iRefBand - (it_stan->poziom); //TEST2
3698 it_stan++;
3699 }
3700}
3701/*****************************************************************************/
3702 void struktura::funkcje_do_pliku_(std::ofstream & plik, char c, double iRefBand, double krok,
3703 double skala) // LUKASZ dodalem cala funkcje
3704{
3705 double szer = prawa.iks - lewa.iks;
3706 double bok = szer / 4;
3707 double x = lewa.iks - bok;
3708 while(x <= lewa.iks)
3709 {
3710 plik << 0.1 * dlugosc_na_A(x) << '\t';
3711 std::vector<stan>::iterator it_stan = rozwiazania.begin();
3712 while(it_stan != rozwiazania.end())
3713 {
3714 if(c == 'e')
3715 plik << (iRefBand + (it_stan->poziom)) +
3716 skala * lewa.funkcjafal(x, it_stan->poziom, it_stan->wspolczynniki[0]);
3717 else if(c == 'h')
3718 plik << (iRefBand - (it_stan->poziom)) +
3719 skala * lewa.funkcjafal(x, it_stan->poziom, it_stan->wspolczynniki[0]);
3720 it_stan++;
3721 if(it_stan != rozwiazania.end())
3722 plik << '\t';
3723 }
3724 plik << '\n';
3725 x += krok;
3726 }
3727 for(int i = 0; i <= (int)kawalki.size() - 1; i++)
3728 {
3729 x = kawalki[i].x_pocz;
3730 while(x <= kawalki[i].x_kon)
3731 {
3732 plik << 0.1 * dlugosc_na_A(x) << '\t';
3733 std::vector<stan>::iterator it_stan = rozwiazania.begin();
3734 while(it_stan != rozwiazania.end())
3735 {
3736 if(c == 'e')
3737 plik << (iRefBand + (it_stan->poziom)) +
3738 skala *
3739 kawalki[i].funkcjafal(x, it_stan->poziom, it_stan->wspolczynniki[2 * i + 1],
3740 it_stan->wspolczynniki[2 * i + 2]);
3741 else if(c == 'h')
3742 plik << (iRefBand - (it_stan->poziom)) +
3743 skala *
3744 kawalki[i].funkcjafal(x, it_stan->poziom, it_stan->wspolczynniki[2 * i + 1],
3745 it_stan->wspolczynniki[2 * i + 2]);
3746 it_stan++;
3747 if(it_stan != rozwiazania.end())
3748 plik << '\t';
3749 }
3750 plik << '\n';
3751 x += krok;
3752 }
3753 }
3754 x = prawa.iks;
3755 while(x <= prawa.iks + bok)
3756 {
3757 plik << 0.1 * dlugosc_na_A(x) << '\t';
3758 std::vector<stan>::iterator it_stan = rozwiazania.begin();
3759 while(it_stan != rozwiazania.end())
3760 {
3761 if(c == 'e')
3762 plik << (iRefBand + (it_stan->poziom)) +
3763 skala * prawa.funkcjafal(x, it_stan->poziom, it_stan->wspolczynniki.back());
3764 else if(c == 'h')
3765 plik << (iRefBand - (it_stan->poziom)) +
3766 skala * prawa.funkcjafal(x, it_stan->poziom, it_stan->wspolczynniki.back());
3767 it_stan++;
3768 if(it_stan != rozwiazania.end())
3769 plik << '\t';
3770 }
3771 plik << '\n';
3772 x += krok;
3773 }
3774}
3775/*****************************************************************************/
3776/*void struktura::funkcje_do_pliku(std::ofstream & plik, double krok)
3777{
3778#ifdef LOGUJ // MD
3779 std::clog << "W f_do_p" << std::endl;
3780#endif
3781 plik << "#\t";
3782 std::vector<stan>::iterator it_stan = rozwiazania.begin();
3783 while(it_stan != rozwiazania.end())
3784 {
3785 plik << " E=" << (it_stan->poziom);
3786 it_stan++;
3787 }
3788 plik << "\n";
3789 double szer = prawa.iks - lewa.iks;
3790 double bok = szer / 4;
3791 double x = lewa.iks - bok;
3792 while(x <= lewa.iks)
3793 {
3794#ifdef MICHAL
3795 plik << dlugosc_na_A(x) << "\t";
3796#else
3797 plik << (0.1 * dlugosc_na_A(x)) << "\t"; // LUKASZ teraz w nm
3798#endif
3799 it_stan = rozwiazania.begin();
3800 while(it_stan != rozwiazania.end())
3801 {
3802 plik << lewa.funkcjafal(x, it_stan->poziom, it_stan->wspolczynniki[0]) << " ";
3803 it_stan++;
3804 }
3805 plik << "\n";
3806 x += krok;
3807 }
3808 for(int i = 0; i <= (int)kawalki.size() - 1; i++)
3809 {
3810 x = kawalki[i].x_pocz;
3811 while(x <= kawalki[i].x_kon)
3812 {
3813#ifdef MICHAL
3814 plik << dlugosc_na_A(x) << "\t";
3815#else
3816 plik << (0.1 * dlugosc_na_A(x)) << "\t"; // LUKASZ teraz w nm
3817#endif
3818 it_stan = rozwiazania.begin();
3819 while(it_stan != rozwiazania.end())
3820 {
3821 plik << kawalki[i].funkcjafal(x, it_stan->poziom, it_stan->wspolczynniki[2 * i + 1],
3822 it_stan->wspolczynniki[2 * i + 2])
3823 << " ";
3824 it_stan++;
3825 }
3826 plik << "\n";
3827 x += krok;
3828 }
3829 }
3830 x = prawa.iks;
3831 while(x <= prawa.iks + bok)
3832 {
3833 plik << dlugosc_na_A(x) << "\t";
3834 it_stan = rozwiazania.begin();
3835 while(it_stan != rozwiazania.end())
3836 {
3837 plik << prawa.funkcjafal(x, it_stan->poziom, it_stan->wspolczynniki.back()) << " ";
3838 it_stan++;
3839 }
3840 plik << "\n";
3841 x += krok;
3842 }
3843}*/
3844/*****************************************************************************/
3845 void struktura::funkcje_do_pliku(std::ofstream & plik, double krok)
3846{
3847 int setw_1 = 13,
3848 setw_2 = 20,
3849 prec_1 = 5,
3850 prec_2 = 10,
3851 ee = 3;
3852
3853 krok = krok / (struktura::przelm * 0.1); // dodalem LUKASZ 2023-09-12
3854
3855 double X = 0.1; // do zamiany A->nm // dodalem X, LUKASZ
3856 std::clog<<"W f_do_p"<<std::endl;
3857 //plik<<"#\t"; // LUKASZ 2023-09-18
3858 plik << setw(setw_1) << "#"; // LUKASZ 2023-09-18
3859 std::vector<stan>::iterator it_stan = rozwiazania.begin();
3860 while(it_stan != rozwiazania.end())
3861 {
3862 plik << std::fixed << std::scientific;
3863 //plik<<" E="<<setw(setw_)<<std::setprecision(precef_)<<(it_stan->poziom); // LUKASZ 2023-09-18
3864 plik << setw(ee) << "E=";
3865 plik << std::fixed << std::scientific;
3866 plik << setw(setw_2-ee) << std::setprecision(prec_2) << (it_stan->poziom); // LUKASZ 2023-09-18
3867 it_stan++;
3868 }
3869 plik<<"\n";
3870 double szer = prawa.iks - lewa.iks;
3871 //double bok = szer/4; // ROBERT 2023-09-18
3872 double bok = szer; // ROBERT 2023-09-18
3873 double x = lewa.iks - bok;
3874 while(x <= lewa.iks)
3875 {
3876 plik << std::fixed << std::scientific;
3877 plik << setw(setw_1) << std::setprecision(prec_1) << dlugosc_na_A(x)*X; // dodalem X, LUKASZ
3878 //plik<<dlugosc_na_A(x)*X/(struktura::przelm*0.1)<<"\t"; // dodalem X, LUKASZ 2023-09-12
3879 it_stan = rozwiazania.begin();
3880 while(it_stan != rozwiazania.end())
3881 {
3882 plik << std::fixed << std::scientific;
3883 plik << setw(setw_2) << std::setprecision(prec_2) << lewa.funkcjafal(x, it_stan->poziom, it_stan->wspolczynniki[0]);//<<" "; // LUKASZ 2023-09-18
3884 it_stan++;
3885 }
3886 plik<<"\n";
3887 x += krok;
3888 }
3889 for(int i = 0; i <= (int) kawalki.size() - 1; i++)
3890 {
3891 x = kawalki[i].x_pocz;
3892 while(x <= kawalki[i].x_kon)
3893 {
3894 plik << setw(setw_1) << std::setprecision(prec_1) << dlugosc_na_A(x)*X; // dodalem X, LUKASZ
3895 //plik<<dlugosc_na_A(x)*X/(struktura::przelm*0.1)<<"\t"; // dodalem X, LUKASZ 2023-09-12
3896 it_stan = rozwiazania.begin();
3897 while(it_stan != rozwiazania.end())
3898 {
3899 plik << std::fixed << std::scientific;
3900 plik << setw(setw_2) << std::setprecision(prec_2) << kawalki[i].funkcjafal(x, it_stan->poziom, it_stan->wspolczynniki[2*i + 1], it_stan->wspolczynniki[2*i + 2]);//<<" "; // LUKASZ 2023-09-18
3901 it_stan++;
3902 }
3903 plik<<"\n";
3904 x += krok;
3905 }
3906 }
3907 x = prawa.iks ;
3908 while(x <= prawa.iks + bok)
3909 {
3910 plik << setw(setw_1) << std::setprecision(prec_1) << dlugosc_na_A(x)*X; // dodalem X, LUKASZ
3911 //plik<<dlugosc_na_A(x)*X/(struktura::przelm*0.1)<<"\t"; // dodalem X, LUKASZ 2023-09-12
3912 it_stan = rozwiazania.begin();
3913 while(it_stan != rozwiazania.end())
3914 {
3915 plik << std::fixed << std::scientific;
3916 plik << setw(setw_2) << std::setprecision(prec_2) << prawa.funkcjafal(x, it_stan->poziom, it_stan->wspolczynniki.back());//<<" "; // LUKASZ 2023-09-18
3917 it_stan++;
3918 }
3919 plik<<"\n";
3920 x += krok;
3921 }
3922}
3923/*****************************************************************************/
3924 void struktura::funkcja1_do_pliku(std::ofstream & plik, stan & st, double krok)
3925{
3926#ifdef LOGUJ // MD
3927 std::clog << "W f_do_p" << std::endl;
3928#endif
3929 plik << "#\t";
3930 plik << " E=" << (st.poziom);
3931 plik << "\n";
3932 double szer = prawa.iks - lewa.iks;
3933 double bok = szer / 4;
3934 double x = lewa.iks - bok;
3935 while(x <= lewa.iks)
3936 {
3937#ifdef MICHAL
3938 plik << dlugosc_na_A(x) << "\t";
3939#else
3940 plik << (0.1 * dlugosc_na_A(x)) << "\t"; // LUKASZ teraz w nm
3941#endif
3942 plik << lewa.funkcjafal(x, st.poziom, st.wspolczynniki[0]) << " ";
3943 plik << "\n";
3944 x += krok;
3945 }
3946 for(int i = 0; i <= (int)kawalki.size() - 1; i++)
3947 {
3948 x = kawalki[i].x_pocz;
3949 while(x <= kawalki[i].x_kon)
3950 {
3951#ifdef MICHAL
3952 plik << dlugosc_na_A(x) << "\t";
3953#else
3954 plik << (0.1 * dlugosc_na_A(x)) << "\t"; // LUKASZ teraz w nm
3955#endif
3956 plik << kawalki[i].funkcjafal(x, st.poziom, st.wspolczynniki[2 * i + 1], st.wspolczynniki[2 * i + 2]) << " ";
3957 plik << "\n";
3958 x += krok;
3959 }
3960 }
3961 x = prawa.iks;
3962 while(x <= prawa.iks + bok)
3963 {
3964#ifdef MICHAL
3965 plik << dlugosc_na_A(x) << "\t";
3966#else
3967 plik << (0.1 * dlugosc_na_A(x)) << "\t"; // LUKASZ teraz w nm
3968#endif
3969 plik << prawa.funkcjafal(x, st.poziom, st.wspolczynniki.back()) << " ";
3970 plik << "\n";
3971 x += krok;
3972 }
3973}
3974/*****************************************************************************/
3975 double struktura::energia_od_k_na_ntym(double k, int nr_war, int n)
3976{
3977 warstwa * war;
3978 if(nr_war == 0)
3979 {
3980 war = &lewa;
3981 }
3982 else
3983 {
3984 if(nr_war == (int)kawalki.size() + 1)
3985 {
3986 war = &kawalki[nr_war - 1];
3987 }
3988 else
3989 {
3990 war = &prawa;
3991 }
3992 }
3993 return war->Eodk(k) + rozwiazania[n].poziom;
3994}
3995/*****************************************************************************
3996 double dE_po_dl(size_t nr, chrop ch)
3997{
3998 double k = sqrt(2*)
3999 double licznik =
4000}
4001*****************************************************************************/
4002// MD
4003 void obszar_aktywny::_obszar_aktywny(struktura * elektron, const std::vector<struktura *> & dziury, double Eg,
4004 const std::vector<double> * DSO, double chropo, double matelem, double Temp)
4005{
4006 int i; // LUKASZ dodalem to
4007 przekr_max = 0.;
4008 pasmo_przew.push_back(elektron);
4009 pasmo_wal = dziury;
4010 chrop = struktura::dlugosc_z_A(chropo);
4011 broad = 0.;
4012 T_ref = Temp;
4013 double dE;
4014 for(/*int*/ i = 0; i <= (int)pasmo_przew.size() - 1; i++) // przesuwa struktury, zeby 0 bylo w lewej
4015 {
4016 dE = -pasmo_przew[i]->lewa.y;
4017 pasmo_przew[i]->przesun_energie(dE);
4018 }
4019 // std::clog<<"\nkonstr aktyw. Po pierwszym for";
4020 for(/*int*/ i = 0; i <= (int)pasmo_wal.size() - 1; i++) // przesuwa struktury, zeby 0 bylo w lewej
4021 {
4022 dE = -pasmo_wal[i]->lewa.y;
4023 pasmo_wal[i]->przesun_energie(dE);
4024 }
4025 // std::clog<<"\nkonstr aktyw. Po drugim for";
4026 Egcc.push_back(0);
4027 Egcv = std::vector<double>(dziury.size(), Eg);
4028 int liczba_war = dziury[0]->kawalki.size() + 2;
4029 if(DSO)
4030 { // MD
4031 DeltaSO = *DSO;
4032 if(DeltaSO.size() != liczba_war)
4033 {
4034 Error err("Niepoprawny rozmiar DeltaSO");
4035 err << ": jest " << DeltaSO.size() << ", powinno byc" << liczba_war;
4036 throw err; // MD
4037 }
4038 }
4039 el_mac.reserve(liczba_war);
4040 for(i = 0; i < liczba_war; i++) // LUKASZ mozna podac wlasny element macierzowy
4041 {
4042 if(matelem <= 0.)
4043 el_mac.push_back(element(i));
4044 else
4045 el_mac.push_back(matelem);
4046#ifdef LOGUJ // MD
4047 std::cerr << "Elem. mac. dla warstwy " << i << ": " << el_mac.back() << "\n"; // LUKASZ info
4048#endif
4049 }
4050 zrob_macierze_przejsc();
4051}
4052/*****************************************************************************/
4053 void obszar_aktywny::_obszar_aktywny(struktura * elektron, const std::vector<struktura *> & dziury,
4054 struktura * elektron_m, const std::vector<struktura *> & dziury_m, double Eg,
4055 const std::vector<double> * DSO, double br, double matelem, double Temp)
4056{
4057 int i; // LUKASZ Visual 6
4058 if(elektron->rozwiazania.size() > elektron_m->rozwiazania.size())
4059 {
4060 Error err;
4061 err << "Za mala Liczba rozwian dla struktury elektronowej w strukturze zmodyfikowanej:\n"
4062 << "zwykle maja " << (elektron->rozwiazania.size()) << " zmodyfikowane maja "
4063 << (elektron_m->rozwiazania.size());
4064 throw err;
4065 }
4066#ifdef LOGUJ // MD
4067 if(elektron->rozwiazania.size() != elektron_m->rozwiazania.size())
4068 std::cerr << "Liczba rozwian dla struktur elektronowych jest rozna!\n"
4069 << "zwykle maja " << (elektron->rozwiazania.size()) << " zmodyfikowane maja "
4070 << (elektron_m->rozwiazania.size());
4071#endif
4072 for(/*int*/ i = 0; i < (int)dziury.size(); ++i) // LUKASZ Visual 6
4073 {
4074 if(dziury[i]->rozwiazania.size() > dziury_m[i]->rozwiazania.size())
4075 {
4076 Error err;
4077 err << "Za mala liczba rozwian dla dziur " << i << " w strukturze zmodyfikowanej:\n"
4078 << "zwykle maja " << (dziury[i]->rozwiazania.size()) << " zmodyfikowane maja "
4079 << (dziury_m[i]->rozwiazania.size());
4080 throw err;
4081 }
4082#ifdef LOGUJ // MD
4083 if(dziury[i]->rozwiazania.size() != dziury_m[i]->rozwiazania.size())
4084 std::cerr << "Liczba rozwian dla dziur " << i << " jest rozna:\n"
4085 << "zwykle maja " << (dziury[i]->rozwiazania.size()) << " zmodyfikowane maja "
4086 << (dziury_m[i]->rozwiazania.size());
4087#endif
4088 }
4089 _obszar_aktywny(elektron, dziury, Eg, DSO, 0., matelem, Temp);
4090 pasmo_przew_mod.push_back(elektron_m);
4091 pasmo_wal_mod = dziury_m;
4092 broad = br;
4093}
4094/*****************************************************************************/
4095// LUKASZ dodalem Dso jako wektor oraz el. mac.
4096 obszar_aktywny::obszar_aktywny(struktura * elektron, const std::vector<struktura *> & dziury, double Eg,
4097 const std::vector<double> & DSO, double chropo, double matelem, double Temp)
4098{
4099 _obszar_aktywny(elektron, dziury, Eg, &DSO, chropo, matelem, Temp);
4100}
4101/*****************************************************************************/
4102// LUKASZ dodalem Dso jako wektor
4103 obszar_aktywny::obszar_aktywny(struktura * elektron, const std::vector<struktura *> & dziury, struktura * elektron_m,
4104 const std::vector<struktura *> dziury_m, double Eg, const std::vector<double> & DSO,
4105 double br, double matelem, double Temp)
4106{
4107 _obszar_aktywny(elektron, dziury, elektron_m, dziury_m, Eg, &DSO, br, matelem, Temp);
4108}
4109/*****************************************************************************/
4110 obszar_aktywny::obszar_aktywny(struktura * elektron, const std::vector<struktura *> & dziury, double Eg, double DSO,
4111 double chropo, double matelem, double Temp)
4112{
4113 DeltaSO.assign(dziury[0]->kawalki.size() + 2, DSO);
4114 _obszar_aktywny(elektron, dziury, Eg, NULL, chropo, matelem, Temp);
4115}
4116/*****************************************************************************/
4117#ifdef MICHAL
4118 obszar_aktywny::obszar_aktywny(struktura * elektron, const std::vector<struktura *> & dziury, struktura * elektron_m,
4119 const std::vector<struktura *> dziury_m, double Eg, double DSO, double br, double Temp,
4120 double matelem) // MD
4121#else
4122 obszar_aktywny::obszar_aktywny(struktura * elektron, const std::vector<struktura *> & dziury, struktura * elektron_m,
4123 const std::vector<struktura *> dziury_m, double Eg, double DSO, double br,
4124 double matelem, double Temp) // MD
4125#endif
4126{
4127 DeltaSO.assign(dziury[0]->kawalki.size() + 2, DSO);
4128 _obszar_aktywny(elektron, dziury, elektron_m, dziury_m, Eg, NULL, br, matelem, Temp);
4129}
4130/*****************************************************************************/
4131 void obszar_aktywny::zapisz_poziomy(std::string nazwa)
4132{
4133 std::ofstream plik;
4134 plik.precision(16);
4135 string nazwa_konk;
4136 for(int ic = 0; ic <= (int)pasmo_przew.size() - 1; ++ic)
4137 {
4138 nazwa_konk = nazwa + "_" + pasmo_przew[ic]->nazwa;
4139 plik.open(nazwa_konk.c_str());
4140 plik << std::scientific;
4141 for(int j = 0; j <= (int)pasmo_przew[ic]->rozwiazania.size() - 1; ++j)
4142 {
4143 plik << pasmo_przew[ic]->rozwiazania[j].poziom << " ";
4144 }
4145 plik << "\n";
4146 plik.close();
4147 if(broad)
4148 {
4149 nazwa_konk = nazwa + "_" + pasmo_przew[ic]->nazwa + "_mod";
4150 plik.open(nazwa_konk.c_str());
4151 plik << std::scientific;
4152 for(int j = 0; j <= (int)pasmo_przew_mod[ic]->rozwiazania.size() - 1; ++j)
4153 {
4154 plik << pasmo_przew_mod[ic]->rozwiazania[j].poziom << " ";
4155 }
4156 plik << "\n";
4157 plik.close();
4158 }
4159 }
4160
4161 for(int iv = 0; iv <= (int)pasmo_wal.size() - 1; ++iv)
4162 {
4163 nazwa_konk = nazwa + "_" + pasmo_wal[iv]->nazwa;
4164 plik.open(nazwa_konk.c_str());
4165 plik << std::scientific;
4166 for(int j = 0; j <= (int)pasmo_wal[iv]->rozwiazania.size() - 1; ++j)
4167 {
4168 plik << pasmo_wal[iv]->rozwiazania[j].poziom << " ";
4169 }
4170 plik << "\n";
4171 plik.close();
4172 if(broad)
4173 {
4174 nazwa_konk = nazwa + "_" + pasmo_wal[iv]->nazwa + "_mod";
4175 plik.open(nazwa_konk.c_str());
4176 plik << std::scientific;
4177 for(int j = 0; j <= (int)pasmo_wal_mod[iv]->rozwiazania.size() - 1; ++j)
4178 {
4179 plik << pasmo_wal_mod[iv]->rozwiazania[j].poziom << " ";
4180 }
4181 plik << "\n";
4182 plik.close();
4183 }
4184 }
4185}
4186/*****************************************************************************/
4187 void obszar_aktywny::paryiprzekrycia_dopliku(std::ofstream & plik, int nr_c, int nr_v)
4188{
4189 struktura * el = pasmo_przew[nr_c];
4190 struktura * dziu = pasmo_wal[nr_v];
4191 A2D * m_prz = calki_przekrycia[nr_c][nr_v];
4192 double E0;
4193 for(int nrpoz_el = 0; nrpoz_el <= int(el->rozwiazania.size()) - 1; nrpoz_el++)
4194 for(int nrpoz_dziu = 0; nrpoz_dziu <= int(dziu->rozwiazania.size()) - 1; nrpoz_dziu++)
4195 {
4196 E0 = Egcv[nr_v] - Egcc[nr_c] + el->rozwiazania[nrpoz_el].poziom + dziu->rozwiazania[nrpoz_dziu].poziom;
4197 plik << E0 << " " << ((*m_prz)[nrpoz_el][nrpoz_dziu]) << "\n";
4198 }
4199}
4200/*****************************************************************************/
4201 double obszar_aktywny::min_przerwa_energetyczna() const
4202{
4203 double przerwa = pasmo_przew[0]->dol + pasmo_wal[0]->dol + Egcv[0];
4204 for(int i = 0; i <= (int)pasmo_przew.size() - 1; i++)
4205 for(int j = 0; j <= (int)pasmo_wal.size() - 1; j++)
4206 {
4207 przerwa = (przerwa > pasmo_przew[i]->dol + pasmo_wal[j]->dol + Egcc[i] + Egcv[j])
4208 ? (pasmo_przew[i]->dol + pasmo_wal[j]->dol + Egcc[i] + Egcv[j])
4209 : przerwa;
4210 }
4211 return przerwa;
4212}
4213/*****************************************************************************/
4214 void obszar_aktywny::policz_calki(const struktura * elektron, const struktura * dziura, A2D & macierz,
4215 TNT::Array2D<std::vector<double> > & wekt_calk_kaw)
4216{
4217 double tymcz;
4218#ifdef LOGUJ // MD
4219 std::cerr << "W funkcji policz_calki. " << elektron->rozwiazania.size() << "x" << dziura->rozwiazania.size() << "\n";
4220#endif
4221 for(int i = 0; i <= (int)elektron->rozwiazania.size() - 1; i++)
4222 for(int j = 0; j <= (int)dziura->rozwiazania.size() - 1; j++)
4223 {
4224 tymcz = calka_ij(elektron, dziura, i, j, wekt_calk_kaw[i][j]);
4225 macierz[i][j] = tymcz * tymcz;
4226 // std::cerr<<"\n w indeksie "<<i<<", "<<j<<" jest "<<macierz[i][j]<<"\n";
4227 if(macierz[i][j] > przekr_max)
4228 {
4229 przekr_max = macierz[i][j];
4230 }
4231 }
4232}
4233/*****************************************************************************
4234 void obszar_aktywny::policz_calki_kawalki(const struktura * elektron, const struktura * dziura,
4235 TNT::Array2D<vector<double> > & macierz)
4236{
4237 double tymcz;
4238 for(int i = 0; i <= (int) elektron->rozwiazania.size() - 1; i++)
4239 for(int j = 0; j <= (int) dziura->rozwiazania.size() - 1; j++)
4240 {
4241 tymcz = calka_ij(elektron, dziura, i, j);
4242 macierz[i][j] = (tymcz < 0)? -tymcz:tymcz;
4243 if(macierz[i][j] > przekr_max)
4244 {
4245 przekr_max = macierz[i][j];
4246 }
4247 }
4248}
4249*****************************************************************************/
4250 double obszar_aktywny::iloczyn_pierwotna_bezpola(double x, int nr_war, const struktura * struk1,
4251 const struktura * struk2, int i, int j)
4252{
4253 double Ec = struk1->rozwiazania[i].poziom;
4254 double Ev = struk2->rozwiazania[j].poziom;
4255 double Ac, Bc, Av, Bv;
4256 double wynik;
4257 if(nr_war == 0) // lewa
4258 {
4259 Bc = struk1->rozwiazania[i].wspolczynniki[0];
4260 Bv = struk2->rozwiazania[j].wspolczynniki[0];
4261 wynik = (struk1->lewa.funkcjafal(x, Ec, Bc) * struk2->lewa.funkcjafal_prim(x, Ev, Bv) -
4262 struk1->lewa.funkcjafal_prim(x, Ec, Bc) * struk2->lewa.funkcjafal(x, Ev, Bv)) /
4263 (struk1->lewa.k_kwadr(Ec) - struk2->lewa.k_kwadr(Ev));
4264 return wynik;
4265 }
4266 if(nr_war == (int)struk1->kawalki.size() + 1) // prawa
4267 {
4268 Ac = struk1->rozwiazania[i].wspolczynniki.back();
4269 Av = struk2->rozwiazania[j].wspolczynniki.back();
4270 wynik = (struk1->prawa.funkcjafal(x, Ec, Ac) * struk2->prawa.funkcjafal_prim(x, Ev, Av) -
4271 struk1->prawa.funkcjafal_prim(x, Ec, Ac) * struk2->prawa.funkcjafal(x, Ev, Av)) /
4272 (struk1->prawa.k_kwadr(Ec) - struk2->prawa.k_kwadr(Ev));
4273 return wynik;
4274 }
4275
4276 Ac = struk1->rozwiazania[i].wspolczynniki[2 * nr_war + 1];
4277 Av = struk2->rozwiazania[j].wspolczynniki[2 * nr_war + 1];
4278 Bc = struk1->rozwiazania[i].wspolczynniki[2 * nr_war + 2];
4279 Bv = struk2->rozwiazania[j].wspolczynniki[2 * nr_war + 2];
4280
4281 wynik = (struk1->kawalki[nr_war].funkcjafal(x, Ec, Ac, Bc) * struk2->kawalki[nr_war].funkcjafal_prim(x, Ev, Av, Bv) -
4282 struk1->kawalki[nr_war].funkcjafal_prim(x, Ec, Ac, Bc) * struk2->kawalki[nr_war].funkcjafal(x, Ev, Av, Bv)) /
4283 (struk1->kawalki[nr_war].k_kwadr(Ec) - struk2->kawalki[nr_war].k_kwadr(Ev));
4284 return wynik;
4285}
4286/*****************************************************************************/
4287 double obszar_aktywny::calka_iloczyn_zpolem(int nr_war, const struktura * struk1, const struktura * struk2, int i,
4288 int j) // numeryczne calkowanie
4289{
4290#ifdef LOGUJ // MD
4291 std::clog << "\nW calk numer. Warstwa " << nr_war << " poziom el " << i << " poziom j " << j << "\n";
4292#endif
4293 double krok = 1.; // na razie krok na sztywno przelm (ok 2.6) A
4294 double Ec = struk1->rozwiazania[i].poziom;
4295 double Ev = struk2->rozwiazania[j].poziom;
4296 double Ac, Bc, Av, Bv;
4297 double wynik = 0;
4298 double x_pocz = struk1->kawalki[nr_war].x_pocz;
4299 double x_kon = struk1->kawalki[nr_war].x_kon;
4300 double szer = x_kon - x_pocz;
4301 int podzial = ceill(szer / krok);
4302 krok = szer / podzial; // wyrownanie kroku
4303 Ac = struk1->rozwiazania[i].wspolczynniki[2 * nr_war + 1];
4304 Av = struk2->rozwiazania[j].wspolczynniki[2 * nr_war + 1];
4305 Bc = struk1->rozwiazania[i].wspolczynniki[2 * nr_war + 2];
4306 Bv = struk2->rozwiazania[j].wspolczynniki[2 * nr_war + 2];
4307 double x = x_pocz + krok / 2;
4308 int ii; // LUKASZ Visual 6
4309 // for(int i = 0; i<= podzial - 1; i++) // LUKASZ dodalem komentarz
4310 for(ii = 0; ii <= podzial - 1; ii++) // LUKASZ dodalem to
4311 {
4312#ifdef LOGUJ // MD
4313 std::clog << "\nwynik = " << wynik << " "; // LUKASZ dodalem komentarz
4314#endif
4315 wynik += struk1->kawalki[nr_war].funkcjafal(x, Ec, Ac, Bc) * struk2->kawalki[nr_war].funkcjafal(x, Ev, Av, Bv);
4316 x += krok;
4317 }
4318 wynik *= krok;
4319 return wynik;
4320}
4321/*****************************************************************************/
4322 double obszar_aktywny::calka_ij(const struktura * elektron, const struktura * dziura, int i, int j,
4323 vector<double> & wektor_calk_kaw)
4324{
4325 double Ec = elektron->rozwiazania[i].poziom;
4326 double Ev = dziura->rozwiazania[j].poziom;
4327 double xk = elektron->lewa.iks;
4328 double Ac, Bc, Av, Bv;
4329 double calk_kaw = 0;
4330 Bc = elektron->rozwiazania[i].wspolczynniki[0];
4331 Bv = dziura->rozwiazania[j].wspolczynniki[0];
4332 double calka = elektron->lewa.funkcjafal(xk, Ec, Bc) * dziura->lewa.funkcjafal_prim(xk, Ev, Bv) -
4333 elektron->lewa.funkcjafal_prim(xk, Ec, Bc) * dziura->lewa.funkcjafal(xk, Ev, Bv);
4334 calka = calka / (elektron->lewa.k_kwadr(Ec) - dziura->lewa.k_kwadr(Ev)); // Taki sprytny wzor na calke mozna dostac
4335 wektor_calk_kaw.push_back(calka);
4336 // std::clog<<" calka w lewej = "<<calka<<"\n";
4337 double pierwk, pierwp, xp;
4338 for(int war = 0; war <= (int)elektron->kawalki.size() - 1; war++)
4339 {
4340 if((elektron->kawalki[war].pole == 0) &&
4341 (dziura->kawalki[war].pole == 0)) // trzeba posprzatac, i wywolywac funkcje tutaj
4342 {
4343 xp = elektron->kawalki[war].x_pocz;
4344 xk = elektron->kawalki[war].x_kon;
4345
4346 Ac = elektron->rozwiazania[i].wspolczynniki[2 * war + 1];
4347 Av = dziura->rozwiazania[j].wspolczynniki[2 * war + 1];
4348 Bc = elektron->rozwiazania[i].wspolczynniki[2 * war + 2];
4349 Bv = dziura->rozwiazania[j].wspolczynniki[2 * war + 2];
4350
4351 pierwk =
4352 elektron->kawalki[war].funkcjafal(xk, Ec, Ac, Bc) * dziura->kawalki[war].funkcjafal_prim(xk, Ev, Av, Bv) -
4353 elektron->kawalki[war].funkcjafal_prim(xk, Ec, Ac, Bc) * dziura->kawalki[war].funkcjafal(xk, Ev, Av, Bv);
4354 pierwp =
4355 elektron->kawalki[war].funkcjafal(xp, Ec, Ac, Bc) * dziura->kawalki[war].funkcjafal_prim(xp, Ev, Av, Bv) -
4356 elektron->kawalki[war].funkcjafal_prim(xp, Ec, Ac, Bc) * dziura->kawalki[war].funkcjafal(xp, Ev, Av, Bv);
4357 calk_kaw = (pierwk - pierwp) / (elektron->kawalki[war].k_kwadr(Ec) - dziura->kawalki[war].k_kwadr(Ev));
4358 wektor_calk_kaw.push_back(calk_kaw);
4359 calka += calk_kaw;
4360 }
4361 else // numerycznie na razie
4362 {
4363 calk_kaw = calka_iloczyn_zpolem(war, elektron, dziura, i, j);
4364 wektor_calk_kaw.push_back(calk_kaw);
4365 calka += calk_kaw;
4366 }
4367 // std::cerr<<"\nwarstwa nr "<<war<<"\tcalka kawalek = "<<calk_kaw;
4368 }
4369 xp = elektron->prawa.iks;
4370
4371 Ac = elektron->rozwiazania[i].wspolczynniki.back();
4372 Av = dziura->rozwiazania[j].wspolczynniki.back();
4373 calk_kaw = -(elektron->prawa.funkcjafal(xp, Ec, Ac) * dziura->prawa.funkcjafal_prim(xp, Ev, Av) -
4374 elektron->prawa.funkcjafal_prim(xp, Ec, Ac) * dziura->prawa.funkcjafal(xp, Ev, Av)) /
4375 (elektron->prawa.k_kwadr(Ec) - dziura->prawa.k_kwadr(Ev)); // -= bo calka jest od xp do +oo
4376 wektor_calk_kaw.push_back(calk_kaw);
4377 calka += calk_kaw;
4378 return calka;
4379}
4380/*****************************************************************************/
4381 void obszar_aktywny::zrob_macierze_przejsc()
4382{
4383#ifdef LOGUJ // MD
4384 std::cerr << "W funkcji zrob_macierze_przejsc\n";
4385#endif
4386 A2D * macierz_calek;
4387 TNT::Array2D<std::vector<double> > * macierz_kawalkow;
4388 calki_przekrycia.resize(pasmo_przew.size());
4389 calki_przekrycia_kawalki.resize(pasmo_przew.size());
4390 for(int i = 0; i <= (int)calki_przekrycia.size() - 1; i++)
4391 {
4392 calki_przekrycia[i].resize(pasmo_wal.size());
4393 calki_przekrycia_kawalki[i].resize(pasmo_wal.size());
4394 }
4395 for(int c = 0; c <= (int)pasmo_przew.size() - 1; c++)
4396 {
4397 for(int v = 0; v <= (int)pasmo_wal.size() - 1; v++)
4398 {
4399 macierz_calek = new A2D(pasmo_przew[c]->rozwiazania.size(), pasmo_wal[v]->rozwiazania.size());
4400 macierz_kawalkow = new TNT::Array2D<std::vector<double> >(pasmo_przew[c]->rozwiazania.size(),
4401 pasmo_wal[v]->rozwiazania.size());
4402 policz_calki(pasmo_przew[c], pasmo_wal[v], *macierz_calek, *macierz_kawalkow);
4403 calki_przekrycia[c][v] = macierz_calek;
4404 calki_przekrycia_kawalki[c][v] = macierz_kawalkow;
4405 // std::cerr<<"Macierz przejsc:"<<(*macierz_calek)<<"\n";
4406 }
4407 }
4408}
4409/*****************************************************************************
4410 void obszar_aktywny::zrob_macierze_kawalkow()
4411{
4412 TNT::Array2D<vector<double> > * macierz_kaw;
4413 calki_przekrycia_kawalki.resize(pasmo_przew.size());
4414 for(int i = 0; i <= (int) calki_przekrycia_kawalki.size() - 1; i++)
4415 {
4416 calki_przekrycia_kawalki[i].resize(pasmo_wal.size());
4417 }
4418 for(int c = 0; c <= (int) pasmo_przew.size() - 1; c++)
4419 {
4420 for(int v = 0; v <= (int) pasmo_wal.size() - 1; v++)
4421 {
4422 macierz = new TNT::Array2D<vector<double> >(pasmo_przew[c]->rozwiazania.size(),
4423 pasmo_wal[v]->rozwiazania.size()); policz_calki_kawalki(pasmo_przew[c], pasmo_wal[v], *macierz);
4424 calki_przekrycia_kawalki[c][v] = macierz;
4425 std::clog<<"Macierz przejsc:"<<(*macierz)<<"\n";
4426 }
4427 }
4428}
4429*****************************************************************************/
4430 double obszar_aktywny::element(int nr_war) // Do przerobienia
4431{
4432 warstwa *warc, *warv;
4433 if(nr_war == 0)
4434 {
4435 warc = &pasmo_przew[0]->lewa;
4436 warv = &pasmo_wal[0]->lewa;
4437 }
4438 else
4439 {
4440 if(nr_war < (int)pasmo_przew[0]->kawalki.size() + 1)
4441 {
4442 warc = &pasmo_przew[0]->kawalki[nr_war - 1];
4443 warv = &pasmo_wal[0]->kawalki[nr_war - 1];
4444 }
4445 else
4446 {
4447 warc = &pasmo_przew[0]->prawa;
4448 warv = &pasmo_wal[0]->prawa;
4449 }
4450 }
4451 double Eg = Egcv[0] + warc->y_pocz + warv->y_pocz;
4452#ifdef LOGUJ // MD
4453 std::cerr << "\nW elemencie: Eg = " << Eg << "\n";
4454#endif
4455 return (1 / warc->masa_p(0.) - 1) * (Eg + DeltaSO[nr_war]) * Eg / (Eg + 2 * DeltaSO[nr_war] / 3) / 2;
4456}
4457/*****************************************************************************/
4458 void obszar_aktywny::ustaw_element(double iM) // dodalem metode LUKASZ LUKI 23.05.2023
4459{
4460 for(int i = 0; i < (int) el_mac.size(); i++)
4461 {
4462 el_mac[i] = iM;
4463 std::cout << el_mac[i] << "\n";
4464 }
4465}
4466/*****************************************************************************/
4467 double wzmocnienie::kodE(double E, double mc, double mv)
4468{
4469 double m = (1 / mc + 1 / mv);
4470 return sqrt(2 * E / m);
4471}
4472/*****************************************************************************/
4473 double wzmocnienie::rored(double, double mc, double mv)
4474{
4475 double m = (1 / mc + 1 / mv);
4476 return 1 / (m * 2 * struktura::pi * szer_do_wzmoc);
4477}
4478/*****************************************************************************/
4479 double wzmocnienie::erf_dorored(double E, double E0, double sigma)
4480{
4481 if(sigma <= 0)
4482 {
4483 Error err; // MD
4484 err << "\nsigma = " << sigma << "!\n";
4485 throw err;
4486 }
4487 return 0.5 * (1 + erf((E - E0) / (sqrt(2) * sigma)));
4488}
4489/*****************************************************************************/
4490 double wzmocnienie::rored_posz(double E, double E0, double mc, double mv,
4491 double sigma) // gestosc do chopowatej studni o nominalnej roznicy energii poziomow E0.
4492 // Wersja najprostsza -- jedno poszerzenie na wszystko
4493{
4494 if(sigma <= 0)
4495 {
4496 Error err; // MD
4497 err << "\nsigma = " << sigma << "!\n";
4498 throw err;
4499 }
4500 double m = (1 / mc + 1 / mv);
4501 // double sigma = posz_en
4502 return erf_dorored(E, E0, sigma) / (m * 2 * struktura::pi * szer_do_wzmoc);
4503}
4504/*****************************************************************************
4505 wzmocnienie::wzmocnienie(obszar_aktywny * obsz, double konc_pow, double temp, double wsp_zal): pasma(obsz)
4506{
4507 // pasma = obsz;
4508 nosniki_c = przel_gest_z_cm2(konc_pow);
4509 nosniki_v = nosniki_c;
4510 T = temp;
4511 ustaw_przerwy(); //ustawia przerwy energetyczne dla danej temperatury
4512 n_r = wsp_zal;
4513 szer_do_wzmoc = pasma->pasmo_przew[0]->kawalki.back().x_kon - pasma->pasmo_przew[0]->kawalki.front().x_pocz;
4514 policz_qFlc();
4515 policz_qFlv();
4516}
4517*****************************************************************************/
4518 wzmocnienie::wzmocnienie(obszar_aktywny * obsz, double konc_pow, double temp, double wsp_zal, double poprawkaEg,
4519 double szdowzm, Wersja wersja)
4520 : pasma(obsz), // szdowzm w A
4521 wersja(wersja) // MD
4522{
4523 // pasma = obsz;
4524 nosniki_c = przel_gest_z_cm2(konc_pow);
4525 nosniki_v = nosniki_c;
4526 T = temp;
4527 ustaw_przerwy(poprawkaEg); // ustawia przerwy energetyczne dla danej temperatury
4528 n_r = wsp_zal;
4529 warstwy_do_nosnikow = std::set<int>();
4530 double szergwiazd = 0.0;
4531 int i, j;
4532 for(i = 0; i <= (int)pasma->pasmo_przew.size() - 1; ++i)
4533 {
4534 for(j = 0; j <= (int)pasma->pasmo_przew[i]->gwiazdki.size() - 1; ++j)
4535 {
4536 warstwy_do_nosnikow.insert(pasma->pasmo_przew[i]->gwiazdki[j]);
4537 }
4538 }
4539 for(i = 0; i <= (int)pasma->pasmo_wal.size() - 1; ++i)
4540 {
4541 for(j = 0; j <= (int)pasma->pasmo_wal[i]->gwiazdki.size() - 1; ++j)
4542 {
4543 warstwy_do_nosnikow.insert(pasma->pasmo_wal[i]->gwiazdki[j]);
4544 }
4545 }
4546 if(warstwy_do_nosnikow.empty()) // nie bylo zadnej gwiazdki, to wszystkie warstwy
4547 {
4548 for(j = 0; j <= (int)pasma->pasmo_przew[0]->kawalki.size() - 1; ++j)
4549 warstwy_do_nosnikow.insert(j);
4550 }
4551 /*
4552 for(std::set<int>::iterator setit = warstwy_do_nosnikow.begin(); setit != warstwy_do_nosnikow.end(); ++setit)
4553 {
4554 std::clog<<"gwiazdka w "<<(*setit)<<"\n";
4555 }
4556 */
4557 for(std::set<int>::iterator setit = warstwy_do_nosnikow.begin(); setit != warstwy_do_nosnikow.end(); ++setit)
4558 {
4559 szergwiazd += pasma->pasmo_przew[0]->kawalki[*setit].x_kon - pasma->pasmo_przew[0]->kawalki[*setit].x_pocz;
4560 // std::clog<<"gwiazdka w "<<(*setit)<<"\n";
4561 }
4562
4563 if(szdowzm < 0)
4564 {
4565 std::cerr<<"\nW ifie\n: krańce = "<<pasma->pasmo_przew[0]->kawalki.back().x_kon<<", "<<pasma->pasmo_przew[0]->kawalki.front().x_pocz<<"\n";
4566 szer_do_wzmoc = pasma->pasmo_przew[0]->kawalki.back().x_kon - pasma->pasmo_przew[0]->kawalki.front().x_pocz;
4567 std::cerr<<"\nW ifie szer_do_wzmoc = "<< szer_do_wzmoc<<"\n";
4568 }
4569 else
4570 szer_do_wzmoc = szdowzm / struktura::przelm;
4571
4572#ifdef LOGUJ // MD
4573 if(abs(szer_do_wzmoc - szergwiazd) / szer_do_wzmoc > 1e-4)
4574 {
4575 std::cerr << "\nszer_do_wzmoc = " << szdowzm << " A\n szergwiad = " << (szergwiazd * struktura::przelm)
4576 << " A!\n";
4577 }
4578 /* else
4579 {
4580 std::clog<<"\n szerokosci gwiazdek i do wzmoc zgadzaja sie z dokladnoscia " <<(abs(szer_do_wzmoc -
4581 szergwiazd)/szer_do_wzmoc)<<"\n";
4582 }
4583 */
4584#endif
4585 policz_qFlc();
4586 policz_qFlv();
4587}
4588/*****************************************************************************/
4589 void wzmocnienie::ustaw_przerwy(double popr)
4590{
4591 Egcv_T.resize(pasma->Egcv.size());
4592 for(size_t i = 0; i <= pasma->Egcv.size() - 1; i++)
4593 {
4594 Egcv_T[i] = pasma->Egcv[i] + popr; // prymitywne przepisanie z obszaru aktywnego
4595 }
4596}
4597/*****************************************************************************/
4598 void wzmocnienie::policz_qFlc()
4599{
4600 // std::clog<<"Szukanie qFlc\n";
4601 double Fp, Fk, krok;
4602 double np, nk;
4603 Fp = -Egcv_T[0]; // 2; // polowa przerwy
4604 Fk = pasma->pasmo_przew[0]->gora; // skrajna bariera
4605 // std::cout << "\tFp: " << Fp << "\n"; // TEST // TEST LUKASZ
4606 // std::cout << "\tFk: " << Fk << "\n"; // TEST // TEST LUKASZ
4607 krok = pasma->pasmo_przew[0]->gora - pasma->pasmo_przew[0]->dol;
4608 // std::cout << "\tkrok: " << krok << "\n"; // TEST // TEST LUKASZ
4609 np = nosniki_w_c(Fp);
4610 nk = nosniki_w_c(Fk);
4611 // std::cout << "\tnp: " << np << "\n"; // TEST // TEST LUKASZ
4612 // std::cout << "\tnk: " << nk << "\n"; // TEST // TEST LUKASZ
4613 if(nosniki_c < np)
4614 {
4615 Error err; // MD
4616 err << "Za malo nosnikow!\n";
4617 throw err;
4618 }
4619 while(nk < nosniki_c)
4620 {
4621 Fp = Fk;
4622 Fk += krok;
4623 nk = nosniki_w_c(Fk);
4624 }
4625 double (wzmocnienie::*fun)(double) = &wzmocnienie::gdzie_qFlc;
4626 // std::clog<<"Sieczne Fermi\n";
4627 qFlc = sieczne(fun, Fp, Fk);
4628}
4629/*****************************************************************************/
4630 void wzmocnienie::policz_qFlv()
4631{
4632 // std::clog<<"Szukanie qFlv\n";
4633 double Fp, Fk, krok;
4634 double np, nk;
4635 Fp = -Egcv_T[0]; // 2; // polowa przerwy
4636 Fk = pasma->pasmo_wal[0]->gora; // skrajna bariera
4637 krok = pasma->pasmo_wal[0]->gora - pasma->pasmo_wal[0]->dol;
4638 np = nosniki_w_v(Fp);
4639 nk = nosniki_w_v(Fk);
4640 if(nosniki_v < np)
4641 {
4642 Error err; // MD
4643 err << "Za malo nosnikow!\n";
4644 throw err;
4645 }
4646 while(nk < nosniki_v)
4647 {
4648 Fp = Fk;
4649 Fk += krok;
4650 nk = nosniki_w_v(Fk);
4651 }
4652 double (wzmocnienie::*fun)(double) = &wzmocnienie::gdzie_qFlv;
4653 // std::clog<<"Sieczne Fermi\n";
4654 qFlv = -sieczne(fun, Fp, Fk); // minus, bo energie sa odwrocone, bo F-D opisuje obsadzenia elektronowe
4655}
4656/*****************************************************************************/
4657 double wzmocnienie::gdzie_qFlc(double E) { return nosniki_w_c(E) - nosniki_c; }
4658/*****************************************************************************/
4659 double wzmocnienie::gdzie_qFlv(double E) { return nosniki_w_v(E) - nosniki_v; }
4660/*****************************************************************************/
4661 double wzmocnienie::nosniki_w_c(double Fl)
4662{
4663 // std::cout << "\tT: " << T << "\n"; // TEST // TEST LUKASZ
4664 double przes;
4665 double n = 0.0;
4666 if(warstwy_do_nosnikow.empty())
4667 n = pasma->pasmo_przew[0]->ilenosnikow(Fl, T);
4668 else
4669 n = pasma->pasmo_przew[0]->ilenosnikow(Fl, T, warstwy_do_nosnikow);
4670 // std::cout << "\tn: " << n << "\n"; // TEST // TEST LUKASZ
4671 for(int i = 1; i <= (int)pasma->pasmo_przew.size() - 1; i++)
4672 {
4673 // std::clog<<"Drugie pasmo c\n";
4674 przes = -pasma->Egcc[i]; // bo dodatnie Egcc oznacza przesuniecie w dol
4675 if(warstwy_do_nosnikow.empty())
4676 n += pasma->pasmo_przew[i]->ilenosnikow(Fl - przes, T);
4677 else
4678 n += pasma->pasmo_przew[i]->ilenosnikow(Fl - przes, T, warstwy_do_nosnikow);
4679 }
4680 // std::clog<<"Fl = "<<Fl<<"\tgest pow w 1/cm^2 = "<<przel_gest_na_cm2(n)<<"\n";
4681 return n;
4682}
4683/*****************************************************************************/
4684 double wzmocnienie::nosniki_w_v(double Fl)
4685{
4686 double przes;
4687 double n = 0.0;
4688 if(warstwy_do_nosnikow.empty())
4689 n = pasma->pasmo_wal[0]->ilenosnikow(Fl, T);
4690 else
4691 n = pasma->pasmo_wal[0]->ilenosnikow(Fl, T, warstwy_do_nosnikow);
4692 for(int i = 1; i <= (int)pasma->pasmo_wal.size() - 1; i++)
4693 {
4694 // std::clog<<"Drugie pasmo v\n";
4695 przes = Egcv_T[0] - Egcv_T[i]; // przerwa Ev_i - Ev_0. dodatnie oznacza ze, v_1 jest blizej c, czyli poziom
4696 // Fermiego trzeba podniesc dla niego (w geometrii struktury)
4697 // std::cerr<<"policzone przes\n";
4698 if(warstwy_do_nosnikow.empty())
4699 n += pasma->pasmo_wal[i]->ilenosnikow(Fl + przes, T);
4700 else
4701 n += pasma->pasmo_wal[i]->ilenosnikow(Fl + przes, T, warstwy_do_nosnikow);
4702 }
4703 // std::clog<<"Fl = "<<Fl<<"\tgest pow w 1/cm^2 = "<<przel_gest_na_cm2(n)<<"\n";
4704 return n;
4705}
4706/*****************************************************************************/
4707 double wzmocnienie::pozFerm_przew() { return qFlc; }
4708/*****************************************************************************/
4709 double wzmocnienie::pozFerm_wal()
4710{
4711 return -qFlv; // w ukladzie, gdzie studnie sa w dol, a bariera na 0.
4712}
4713/*****************************************************************************/
4714 double wzmocnienie::rozn_poz_Ferm() { return Egcv_T[0] + qFlc - qFlv; }
4715/*****************************************************************************/
4716 double wzmocnienie::szerdowzmoc()
4717{
4718 return szer_do_wzmoc * struktura::przelm; // przelicza na A
4719}
4720/*****************************************************************************/
4721std::vector<double> wzmocnienie::koncentracje_elektronow_w_warstwach() // zwraca w 1/cm^3
4722{
4723 std::vector<double> konce, koncetym;
4724 std::vector<struktura *>::const_iterator piter = pasma->pasmo_przew.begin();
4725 konce = (*piter)->koncentracje_w_warstwach(qFlc, T);
4726 ++piter;
4727 while(piter != pasma->pasmo_przew.end())
4728 {
4729 koncetym = (*piter)->koncentracje_w_warstwach(qFlc, T);
4730 ++piter;
4731 for(int i = 0; i <= (int)konce.size() - 1; ++i)
4732 {
4733 konce[i] += koncetym[i];
4734 }
4735 }
4736 for(int i = 0; i <= (int)konce.size() - 1; ++i)
4737 {
4738 konce[i] = struktura::koncentracja_na_cm_3(konce[i]);
4739 }
4740 return konce;
4741}
4742/*****************************************************************************/
4743std::vector<double> wzmocnienie::koncentracje_dziur_w_warstwach()
4744{
4745 std::vector<double> konce, koncetym;
4746 std::vector<struktura *>::const_iterator piter = pasma->pasmo_wal.begin();
4747 konce =
4748 (*piter)->koncentracje_w_warstwach(-qFlv, T); // minus, bo struktury nie sa swiadome obrocenia energii dla dziur
4749 ++piter;
4750 while(piter != pasma->pasmo_wal.end())
4751 {
4752 koncetym = (*piter)->koncentracje_w_warstwach(-qFlv, T);
4753 ++piter;
4754 for(int i = 0; i <= (int)konce.size() - 1; ++i)
4755 {
4756 konce[i] += koncetym[i];
4757 }
4758 }
4759 for(int i = 0; i <= (int)konce.size() - 1; ++i)
4760 {
4761 konce[i] = struktura::koncentracja_na_cm_3(konce[i]);
4762 }
4763 return konce;
4764}
4765
4766/* double tylenosnikow = 0;
4767 double niepomnozone;
4768 double calkazFD, calkaFD12;
4769 double szer;
4770 double gam32 = sqrt(struktura::pi)/2; // Gamma(3/2)
4771 double kT = kB*T;
4772 double Egorna = lewa.y;
4773 calkaFD12 = gsl_sf_fermi_dirac_half ((qFlc-Egorna)/kT);
4774 std::vector<double> koncentr(kawalki.size() + 2, 0.0);
4775 std::vector<stan>::iterator it;
4776 std::vector<struktura *> pasprz = pasma.pasmo_przew;
4777 // std::clog<<"Liczba stanow = "<<rozwiazania.size()<<"\n";
4778 std::vector<struktura *>::iterator piter;
4779 for(piter = pasprz.begin(); piter != pasprz.begin(); ++piter)
4780 {
4781 koncentr[0] += sqrt(2*piter->lewa.masa_p)*piter->lewa.masa_r *
4782 kT*gam32*sqrt(kT)*2*calkaFD12/(2*struktura::pi*struktura::pi); for(int i = 0; i<= (int) piter->kawalki.size() - 1; ++i)
4783// Suma po warstwach
4784 {
4785 it = piter->rozwiazania.end();
4786 std::clog<<"i = "<<i<<"\n";
4787 tylenosnikow = 0;
4788 while(it != piter->rozwiazania.begin()) // suma po poziomach
4789 {
4790 --it;
4791 if(i == 1)
4792 {
4793 std::clog<<"Zawartosc dla poziomu "<<(it->poziom)<<" wynosi
4794"<<(piter->kawalki[i].norma_kwadr(it->poziom, it->wspolczynniki[2*i + 1], it->wspolczynniki[2*i + 2]))<<"\n";
4795 }
4796 calkazFD = kT*log(1+exp((qFl - it->poziom)/kT));
4797 niepomnozone = piter->kawalki[i].norma_kwadr(it->poziom, it->wspolczynniki[2*i + 1], it->wspolczynniki[2*i
4798+ 2]) * piter->kawalki[i].masa_r; tylenosnikow += niepomnozone*calkazFD/struktura::pi; // spiny juz sa
4799 }
4800 szer = piter->kawalki[i].x_kon - piter->kawalki[i].x_pocz;
4801 koncentr[i + 1] += tylenosnikow/szer + sqrt(2*piter->kawalki[i].masa_p(Egorna))*piter->kawalki[i].masa_r*
4802 kT*gam32*sqrt(kT)*2*calkaFD12/(2*struktura::pi*struktura::pi);
4803 }
4804 }
4805 koncentr.back() = koncentr.front();
4806 return koncentr;
4807}
4808*****************************************************************************/
4809 double wzmocnienie::sieczne(double (wzmocnienie::*f)(double), double pocz, double kon)
4810{
4811#ifdef LOGUJ // MD
4812 std::clog.precision(12);
4813#endif
4814 // const double eps = 1e-14; // limit zmian x
4815 double dokl = 1e-6;
4816 double xp = kon;
4817 double xl = pocz;
4818 /*
4819 if((this->*f)(pocz)*(this->*f)(kon) > 0)
4820 {
4821 Error err; // MD
4822 err<<"Zle krance przedzialu!\n";
4823 throw err;
4824 }
4825 */
4826 double x, fc, fp, fl, xlp, xpp; // -p -- poprzednie krance
4827 fl = (this->*f)(xl);
4828 fp = (this->*f)(xp);
4829 xlp = (xl + xp) / 2;
4830 xpp = xlp; // zeby na pewno bylo na poczatku rozne od prawdzwych koncow
4831 do
4832 {
4833 x = xp - fp * (xp - xl) / (fp - fl);
4834 fc = (this->*f)(x);
4835 if(fc == 0)
4836 {
4837 break;
4838 }
4839 if(fc * fl < 0) // c bedzie nowym prawym koncem
4840 {
4841 // std::clog<<"xlp - xl = "<<(xlp - xl)<<"\n";
4842 if(xlp == xl) // trzeba poprawic ten kraniec (metoda Illinois)
4843 {
4844 // std::clog<<"Lewy Illinois\n";
4845 fl = fl / 2;
4846 }
4847 xpp = xp;
4848 xp = x;
4849 xlp = xl;
4850 fp = fc;
4851 }
4852 else // c bedzie nowym lewym koncem
4853 {
4854 // std::clog<<"xpp - xp = "<<(xpp - xp)<<"\n";
4855 if(xpp == xp) // trzeba poprawic ten kraniec (metoda Illinois)
4856 {
4857 // std::clog<<"Prawy Illinois\n";
4858 fp = fp / 2;
4859 }
4860 xlp = xl;
4861 xl = x;
4862 xpp = xp;
4863 fl = fc;
4864 // fp = (this->*f)(kon);
4865 }
4866 // std::clog<<"x = "<<x<<"\tf(x) = "<<fc<<"\txl = "<<xl<<" xp = "<<xp<<"\n";//<<"\txp - x = "<<(xp -
4867 // x)<<"\n";
4868 }
4869 while(xp - xl >= dokl);
4870 return x;
4871}
4872/*****************************************************************************/
4873 double wzmocnienie::L(double x, double b) { return 1 / (M_PI * b) / (1 + x / b * x / b); }
4874/*****************************************************************************/
4875 double wzmocnienie::wzmocnienie_calk_ze_splotem(double E, double b, double polar,
4876 double blad) // podzial na kawalek o promieniu Rb wokol 0 i reszte
4877{
4878 // double blad = 0.005;
4879 // b energia do poszerzenia w lorentzu
4880 struktura * el = pasma->pasmo_przew[0];
4881 struktura * dziu = pasma->pasmo_wal[0];
4882 double E0pop = Egcv_T[0] - pasma->Egcc[0] + el->rozwiazania[0].poziom +
4883 dziu->rozwiazania[0].poziom; // energia potencjalna + energia prostopadla
4884 double E0min = E0pop;
4885 ;
4886 for(int nr_c = 0; nr_c <= (int)pasma->pasmo_przew.size() - 1; nr_c++)
4887 {
4888 el = pasma->pasmo_przew[nr_c];
4889 for(int nr_v = 0; nr_v <= (int)pasma->pasmo_wal.size() - 1; nr_v++)
4890 {
4891 dziu = pasma->pasmo_wal[nr_v];
4892 // MD E0min = Egcv_T[nr_c] - pasma->Egcc[nr_v] + el->rozwiazania[0].poziom + dziu->rozwiazania[0].poziom;
4893 E0min = Egcv_T[nr_v] - pasma->Egcc[nr_c] + el->rozwiazania[0].poziom + dziu->rozwiazania[0].poziom;
4894 E0min = (E0pop >= E0min) ? E0min : E0pop;
4895 }
4896 }
4897 double a = 2 * (E0min - pasma->min_przerwa_energetyczna()) * pasma->chrop;
4898 // maksima (oszacowne z gory) kolejnych pochodnych erfc(x/a)
4899 double em = 2.;
4900 double epm = 1.13 / a;
4901 double eppm = 1. / (a * a);
4902 double epppm = 2.5 / (a * a * a);
4903 double eppppm = 5 / (a * a * a * a);
4904 // maxima (oszacowne z gory) kolejnych pochodnych lorentza (x/b), pomnozeone przez b
4905 double lm = 1 / M_PI;
4906 double lpm = 0.2 / b;
4907 double lppm = 0.7 / (b * b);
4908 double lpppm = 1.5 / (b * b * b);
4909 double lppppm = 24 / M_PI / (b * b * b * b);
4910 double czwpoch0b = (em * lppppm + 4 * epm * lpppm + 6 * eppm * lppm + 4 * epppm * lpm +
4911 eppppm * lm); // szacowanie (grube) czwartej pochodnej dla |x| < 1, pomnozone przez b
4912 double R = 3.; // w ilu b jest zmiana zageszczenia
4913 double jedenplusR2 = 1 + R * R;
4914 double lmR =
4915 1 / (M_PI * jedenplusR2); // oszacowania modulow kolejnych pochodnych przez wartosc w R, bo dla R >=2 sa malejace
4916 double lpmR = 2 * R / (M_PI * jedenplusR2 * jedenplusR2) / b;
4917 double lppmR = (6 * R * R - 2) / (M_PI * jedenplusR2 * jedenplusR2 * jedenplusR2) / (b * b);
4918 double lpppmR = 24 * R * (R * R - 1) / (M_PI * jedenplusR2 * jedenplusR2 * jedenplusR2 * jedenplusR2) / (b * b * b);
4919 double lppppmR = (120 * R * R * R * R - 240 * R * R + 24) /
4920 (M_PI * jedenplusR2 * jedenplusR2 * jedenplusR2 * jedenplusR2 * jedenplusR2) / (b * b * b * b);
4921 double czwpochRb = (em * lppppmR + 4 * epm * lpppmR + 6 * eppm * lppmR + 4 * epppm * lpmR + eppppm * lmR);
4922
4923 int n0 = pow(2 * R, 5. / 4) * b * pow(czwpoch0b / (180 * blad), 0.25);
4924 int nR = pow((32 - R), 5. / 4) * b * pow(czwpochRb / (180 * blad), 0.25);
4925 if(n0 % 2) // ny powinny byc parzyste
4926 {
4927 n0 += 1;
4928 }
4929 else
4930 {
4931 n0 += 2;
4932 }
4933 if(nR % 2) // ny powinny byc parzyste
4934 {
4935 nR += 1;
4936 }
4937 else
4938 {
4939 nR += 2;
4940 }
4941 // nR *= 2; // chwilowo, do testow
4942 double szer = 2 * R * b;
4943 double h = szer / n0;
4944 double x2j, x2jm1, x2jm2;
4945 double calka1 = 0.;
4946 int j;
4947 for(j = 1; j <= n0 / 2; j++) // w promieniu Rb
4948 {
4949 x2j = -R * b + 2 * j * h;
4950 x2jm1 = x2j - h;
4951 x2jm2 = x2jm1 - h;
4952 calka1 += L(x2jm2, b) * wzmocnienie_calk_bez_splotu(E - x2jm2, polar) +
4953 4 * L(x2jm1, b) * wzmocnienie_calk_bez_splotu(E - x2jm1, polar) + L(x2j, b) * wzmocnienie_calk_bez_splotu(E - x2j, polar);
4954 }
4955 calka1 *= h / 3;
4956 // dla -32b < x -Rb
4957 szer = (32 - R) * b;
4958 h = szer / nR;
4959 double calka2 = 0.;
4960 for(j = 1; j <= nR / 2; j++) // ujemne pol
4961 {
4962 x2j = -32 * b + 2 * j * h;
4963 x2jm1 = x2j - h;
4964 x2jm2 = x2jm1 - h;
4965 calka2 += L(x2jm2, b) * wzmocnienie_calk_bez_splotu(E - x2jm2, polar) +
4966 4 * L(x2jm1, b) * wzmocnienie_calk_bez_splotu(E - x2jm1, polar) + L(x2j, b) * wzmocnienie_calk_bez_splotu(E - x2j, polar);
4967 }
4968 for(j = 1; j <= nR / 2; j++) // dodatnie pol
4969 {
4970 x2j = R * b + 2 * j * h;
4971 x2jm1 = x2j - h;
4972 x2jm2 = x2jm1 - h;
4973 calka2 += L(x2jm2, b) * wzmocnienie_calk_bez_splotu(E - x2jm2, polar) +
4974 4 * L(x2jm1, b) * wzmocnienie_calk_bez_splotu(E - x2jm1, polar) + L(x2j, b) * wzmocnienie_calk_bez_splotu(E - x2j, polar);
4975 }
4976 calka2 *= h / 3;
4977 double calka = calka1 + calka2;
4978#ifdef LOGUJ // MD
4979 std::clog << "\na = " << a << "\t4poch = " << czwpoch0b << "\tn0 = " << n0 << "\tnR = " << nR << "\tcalka = " << calka
4980 << "\n";
4981#endif
4982 return calka;
4983}
4984/*****************************************************************************/
4985 double wzmocnienie::wzmocnienie_od_pary_pasm(double E, size_t nr_c, size_t nr_v, double polar)
4986{
4987 // std::cerr<<"\npasmo walencyjna nr "<<nr_v<<"\n";
4988 struktura * el = pasma->pasmo_przew[nr_c];
4989 struktura * dziu = pasma->pasmo_wal[nr_v];
4990 A2D * m_prz = pasma->calki_przekrycia[nr_c][nr_v];
4991 double wzmoc = 0;
4992 // double minimalna_przerwa = Egcv_T[nr_v] - pasma->Egcc[nr_c] + el->dol + dziu->dol;
4993 // double min_E0 = Egcv_T[nr_v] - pasma->Egcc[nr_c] + el->rozwiazania[0].poziom + dziu->rozwiazania[0].poziom;
4994 // double posz_en = 2*(min_E0 - minimalna_przerwa)*pasma->chrop;
4995 double E0;
4996 double posz_en;
4997 for(int nrpoz_el = 0; nrpoz_el <= int(el->rozwiazania.size()) - 1; nrpoz_el++)
4998 for(int nrpoz_dziu = 0; nrpoz_dziu <= int(dziu->rozwiazania.size()) - 1; nrpoz_dziu++)
4999 {
5000 // std::cerr<<"\nprzekrycie w "<<nrpoz_el<<", "<<nrpoz_dziu
5001 // <<" = "<<(*m_prz)[nrpoz_el][nrpoz_dziu];
5002 // std::cout << "kkk " << Egcv_T[nr_v] << " " << pasma->Egcc[nr_c] << " " << el->rozwiazania[nrpoz_el].poziom <<
5003 // " " << dziu->rozwiazania[nrpoz_dziu].poziom << "\n";
5004 E0 = Egcv_T[nr_v] - pasma->Egcc[nr_c] + el->rozwiazania[nrpoz_el].poziom + dziu->rozwiazania[nrpoz_dziu].poziom;
5005 // MD switch-case (lepsze niż if)
5006 switch(wersja)
5007 {
5008 case Z_CHROPOWATOSCIA:
5009 posz_en = posz_z_chrop(nr_c, nrpoz_el, nr_v, nrpoz_dziu);
5010 break;
5011 case Z_POSZERZENIEM:
5012 posz_en = posz_z_br(nr_c, nrpoz_el, nr_v, nrpoz_dziu);
5013 break;
5014 }
5015 // std::cout << "\tgain from pair of bands (if) - " << ">0?: " << ((*m_prz)[nrpoz_el][nrpoz_dziu] - 0.005) << ",
5016 // >0?: " << (E-E0 + 8*posz_en) << ", nr_c: " << nr_c << ", nr_v: " << nr_v << ", nrpoz_el: " << nrpoz_el << ",
5017 // nrpoz_dziu: " << nrpoz_dziu << "\n"; // LUKASZ dodalem linie
5018 if(((*m_prz)[nrpoz_el][nrpoz_dziu] > 0.005) && (E - E0 > -8 * posz_en)) // czy warto tracic czas
5019 {
5020 // std::cout << "\t\tcalc gain from pair of levels (if)\n"; // LUKASZ dodalem linie
5021 wzmoc += wzmocnienie_od_pary_poziomow(E, nr_c, nrpoz_el, nr_v, nrpoz_dziu, polar);
5022 // std::cerr<<"\nnowy wzmoc = "<<wzmoc;
5023 }
5024 }
5025 return wzmoc;
5026}
5027/*****************************************************************************/
5028 double wzmocnienie::spont_od_pary_pasm(double E, size_t nr_c, size_t nr_v,
5029 double polar) // polar = 0 to TE, polar 1 to TM
5030{
5031 // std::cerr<<"\npasmo walencyjna nr "<<nr_v<<"\n";
5032 struktura * el = pasma->pasmo_przew[nr_c];
5033 struktura * dziu = pasma->pasmo_wal[nr_v];
5034 A2D * m_prz = pasma->calki_przekrycia[nr_c][nr_v];
5035 double spont = 0;
5036 // double minimalna_przerwa = Egcv_T[nr_v] - pasma->Egcc[nr_c] + el->dol + dziu->dol;
5037 // double min_E0 = Egcv_T[nr_v] - pasma->Egcc[nr_c] + el->rozwiazania[0].poziom + dziu->rozwiazania[0].poziom;
5038 double posz_en;
5039 double E0;
5040 for(int nrpoz_el = 0; nrpoz_el <= int(el->rozwiazania.size()) - 1; nrpoz_el++)
5041 for(int nrpoz_dziu = 0; nrpoz_dziu <= int(dziu->rozwiazania.size()) - 1; nrpoz_dziu++)
5042 {
5043 // std::cerr<<"\nprzekrycie w "<<nrpoz_el<<", "<<nrpoz_dziu
5044 // <<" = "<<(*m_prz)[nrpoz_el][nrpoz_dziu];
5045 E0 = Egcv_T[nr_v] - pasma->Egcc[nr_c] + el->rozwiazania[nrpoz_el].poziom + dziu->rozwiazania[nrpoz_dziu].poziom;
5046 // MD switch-case (lepsze niż if)
5047 switch(wersja)
5048 {
5049 case Z_CHROPOWATOSCIA:
5050 posz_en = posz_z_chrop(nr_c, nrpoz_el, nr_v, nrpoz_dziu);
5051 break;
5052 case Z_POSZERZENIEM:
5053 posz_en = posz_z_br(nr_c, nrpoz_el, nr_v, nrpoz_dziu);
5054 break;
5055 }
5056 // std::cerr<<"\n"<<nrpoz_el<<" "<<nrpoz_dziu<<" posz_en = "<<posz_en;
5057 // std::cout << "\tlumi from pair of bands (if) - " << ">0?: " << ((*m_prz)[nrpoz_el][nrpoz_dziu] - 0.005) <<
5058 // ",>0?: " << (E-E0 + 8*posz_en) << ", nr_c: " << nr_c << ", nr_v: " << nr_v << ", nrpoz_el: " << nrpoz_el <<
5059 // ",nrpoz_dziu: " << nrpoz_dziu << "\n"; // LUKASZ dodalem linie
5060 if(((*m_prz)[nrpoz_el][nrpoz_dziu] > 0.005) && (E - E0 > -8 * posz_en)) // czy warto tracic czas
5061 {
5062 // std::cout << "\t\tcalc lumi from pair of levels (if)\n"; // LUKASZ dodalem linie
5063 spont += spont_od_pary_poziomow(E, nr_c, nrpoz_el, nr_v, nrpoz_dziu, polar);
5064 // std::cerr<<"\nnowy wzmoc = "<<wzmoc;
5065 }
5066 }
5067 return spont;
5068}
5069/*****************************************************************************/
5070 double wzmocnienie::wzmocnienie_od_pary_poziomow(double E, size_t nr_c, int poz_c, size_t nr_v, int poz_v, double polar)
5071{
5072 // std::cout << "\t\tgain from pair of levels (0) - " << "nr_c: " << nr_c << ", nr_v: " << nr_v << ", poz_c: " <<
5073 // poz_c << ", poz_v: " << poz_v << "\n";
5074 double wynik;
5075 double cos2tet; // zmiana elementu macierzowego z k_prost
5076 struktura * el = pasma->pasmo_przew[nr_c];
5077 struktura * dziu = pasma->pasmo_wal[nr_v];
5078 double Eg; // lokalna przerwa energetyczna
5079 double E0 = Egcv_T[nr_v] - pasma->Egcc[nr_c] + el->rozwiazania[poz_c].poziom +
5080 dziu->rozwiazania[poz_v].poziom; // energia potencjalna + energia prostopadla
5081 // std::cerr<<"\npoziom_el = "<<el->rozwiazania[poz_c].poziom<<"\n"
5082 // <<"\n\nE = "<<E<<" poziom c = "<<poz_c<<" poziom v = "<<poz_v<<" E0 = "<<E0<<"\n";
5083 double przekr_w_war;
5084 // std::cerr<<"\nTyp dziury = "<<dziu->typ;
5085 std::vector<double> calki_kawalki;
5086
5087 // Usrednianie masy efektywnej
5088 // std::cerr<<"\n prawd = "<<el->rozwiazania[poz_c].prawdopodobienstwa[0];
5089 double srednia_masa_el = el->rozwiazania[poz_c].prawdopodobienstwa[0] * el->lewa.masa_r;
5090 double srednia_masa_dziu = dziu->rozwiazania[poz_v].prawdopodobienstwa[0] * dziu->lewa.masa_r;
5091 int i;
5092 for(i = 0; i <= (int)el->kawalki.size() - 1; i++)
5093 {
5094 // std::cerr<<"\n prawd = "<<el->rozwiazania[poz_c].prawdopodobienstwa[i+1];
5095 srednia_masa_el += el->rozwiazania[poz_c].prawdopodobienstwa[i + 1] * el->kawalki[i].masa_r;
5096 srednia_masa_dziu += dziu->rozwiazania[poz_v].prawdopodobienstwa[i + 1] * dziu->kawalki[i].masa_r;
5097 }
5098 int ost_ind = el->kawalki.size() + 1;
5099 // std::cerr<<"\n prawd = "<<el->rozwiazania[poz_c].prawdopodobienstwa[ost_ind];
5100 srednia_masa_el += el->rozwiazania[poz_c].prawdopodobienstwa[ost_ind] * el->prawa.masa_r;
5101 srednia_masa_dziu += dziu->rozwiazania[poz_v].prawdopodobienstwa[ost_ind] * dziu->prawa.masa_r;
5102 // std::cerr<<"\nSrednie masy:\n elektron = "<< srednia_masa_el<<"\ndziura = "<<srednia_masa_dziu<<"\n";;
5103 double mnoznik_pol;
5104 // koniec usredniania masy
5105
5106 // double srednie_k = (E-E0>0)?kodE(E-E0, srednia_masa_el, srednia_masa_dziu):0.;
5107 double srednie_k_zeznakiem =
5108 (E - E0 > 0) ? kodE(E - E0, srednia_masa_el, srednia_masa_dziu) : -kodE(E0 - E, srednia_masa_el, srednia_masa_dziu);
5109
5110 // double minimalna_przerwa = Egcv_T[nr_v] - pasma->Egcc[nr_c] + el->dol + dziu->dol;
5111 // double min_E0 = Egcv_T[nr_v] - pasma->Egcc[nr_c] + el->rozwiazania[0].poziom + dziu->rozwiazania[0].poziom;
5112 // double posz_en = 2*(min_E0 - minimalna_przerwa)*pasma->chrop; // oszacowanie rozmycia poziomow z powodu
5113 // chropowatosci
5114 double posz_en; // = posz_z_br(nr_c, poz_c, nr_v, poz_v); // LUKASZ dodalem komentarz
5115 // MD switch-case (lepsze niż if)
5116 switch(wersja)
5117 {
5118 case Z_CHROPOWATOSCIA:
5119 posz_en = posz_z_chrop(nr_c, poz_c, nr_v, poz_v);
5120 break;
5121 case Z_POSZERZENIEM:
5122 posz_en = posz_z_br(nr_c, poz_c, nr_v, poz_v);
5123 break;
5124 }
5125
5126 double sr_E_E0_dod = posz_en / (sqrt(2 * struktura::pi)) * exp(-(E - E0) * (E - E0) / (2 * posz_en * posz_en)) +
5127 (E - E0) * erf_dorored(E, E0, posz_en); // srednia energia kinetyczna w plaszczyznie
5128 Eg = Egcv_T[nr_v] - pasma->Egcc[nr_c];
5129 // std::cerr<<"lewa Eg = "<<Eg<<"\n";
5130 // cos2tet= (E0>Eg && E > E0)?(E0-Eg)/(E-Eg):1.0;
5131 cos2tet = (E0 > Eg) ? (E0 - Eg) / (sr_E_E0_dod + E0 - Eg) : 1.;
5132 // std::cerr<<"\ncos2tet = "<<cos2tet<<"\n";
5133 calki_kawalki = (*(pasma->calki_przekrycia_kawalki[nr_c][nr_v]))[poz_c][poz_v];
5134 // std::cerr<<"lewa po calkikawalki\n";
5135 przekr_w_war = calki_kawalki[0];
5136 // std::cerr<<"lewa przed wynikiem\n";
5137 mnoznik_pol = (dziu->typ == struktura::hh) ? (1 + cos2tet + polar * (1 - 3 * cos2tet)) / 2
5138 : (5 - 3 * cos2tet + 3 * polar * (3 - cos2tet)) / 6;
5139 // if (dziu->typ == struktura::hh) std::cout << "\t\tgain from pair of levels (1) - struktura: hh\n"; // LUKASZ
5140 // else std::cout << "\t\tgain from pair of levels (1) - struktura: lh\n"; // LUKASZ dodalem linie
5141 wynik = sqrt(pasma->el_mac[0] * mnoznik_pol) * przekr_w_war;
5142 // std::cerr<<"\nprzekr_w_war = "<<przekr_w_war<<" el_mac = "<<pasma->el_mac[0]<<" wynik = "<<wynik;
5143 for(i = 0; i <= (int)el->kawalki.size() - 1; i++)
5144 {
5145 Eg = Egcv_T[nr_v] - pasma->Egcc[nr_c] + el->kawalki[i].y_pocz +
5146 dziu->kawalki[i].y_pocz; // y_pocz na szybko, moze co innego powinno byc
5147 // cos2tet= (E0>Eg && E > E0)?(E0-Eg)/(E-Eg):1.0;
5148 cos2tet = (E0 > Eg) ? (E0 - Eg) / (sr_E_E0_dod + E0 - Eg) : 1.;
5149 mnoznik_pol = (dziu->typ == struktura::hh) ? (1 + cos2tet + polar * (1 - 3 * cos2tet)) / 2
5150 : (5 - 3 * cos2tet + 3 * polar * (3 - cos2tet)) / 6;
5151 // if (dziu->typ == struktura::hh) std::cout << "\t\tgain from pair of levels (for) - struktura: hh\n"; // LUKASZ
5152 // else std::cout << "\t\tgain from pair of levels (for) - struktura: lh\n"; // LUKASZ dodalem linie
5153 // std::cerr<<"\nkawalek "<<i
5154 // <<" mnoz_pol = "<<mnoznik_pol<<" cos2tet = "<<cos2tet;
5155
5156 przekr_w_war = calki_kawalki[i + 1];
5157 wynik += sqrt(pasma->el_mac[i + 1] * mnoznik_pol) * przekr_w_war;
5158 // std::cerr<<" przekr_w_war = "<<przekr_w_war<< " wynik = "<<wynik;
5159 // std::cerr<<"\nprzekr_w_war = "<<przekr_w_war<<" el_mac = "<<pasma->el_mac[i+1]<<" wynik = "<<wynik;
5160 }
5161 przekr_w_war = calki_kawalki.back();
5162 double energia_elektronu =
5163 el->rozwiazania[poz_c].poziom + srednie_k_zeznakiem * abs(srednie_k_zeznakiem) / (2 * srednia_masa_el);
5164 double energia_dziury =
5165 dziu->rozwiazania[poz_v].poziom + srednie_k_zeznakiem * abs(srednie_k_zeznakiem) / (2 * srednia_masa_dziu);
5166 double rozn_obsadzen = fc(energia_elektronu - pasma->Egcc[nr_c]) - fv(-energia_dziury + Egcv_T[0] - Egcv_T[nr_v]);
5167 // std::cerr<<"\nen_el = "<<energia_elektronu<<" en_dziu = "<<energia_dziury<<" rozn_obsadzen =
5168 // "<<rozn_obsadzen<<"\n";
5169 Eg = Egcv_T[nr_v] - pasma->Egcc[nr_c];
5170 // cos2tet= (E0>Eg && E > E0)?(E0-Eg)/(E-Eg):1.0;
5171 cos2tet = (E0 > Eg) ? (E0 - Eg) / (sr_E_E0_dod + E0 - Eg) : 1.;
5172 mnoznik_pol = (dziu->typ == struktura::hh) ? (1 + cos2tet + polar * (1 - 3 * cos2tet)) / 2
5173 : (5 - 3 * cos2tet + 3 * polar * (3 - cos2tet)) / 6;
5174 // if (dziu->typ == struktura::hh) std::cout << "\t\tgain from pair of levels (2) - struktura: hh\n"; // LUKASZ
5175 // else std::cout << "\t\tgain from pair of levels (2) - struktura: lh\n"; // LUKASZ dodalem linie
5176 wynik += sqrt(pasma->el_mac.back() * mnoznik_pol) * przekr_w_war;
5177 // std::cerr<<"\nprzekr_w_war = "<<przekr_w_war<<" el_mac = "<<pasma->el_mac.back()<<" wynik = "<<wynik<<" rored =
5178 // "<<rored(srednie_k, srednia_masa_el, srednia_masa_dziu)<<"\n";
5179 wynik *= wynik; // dopiero teraz kwadrat modulu
5180 // std::cerr<<"\nwynik = "<<wynik;
5181
5182 wynik *= rored_posz(E, E0, srednia_masa_el, srednia_masa_dziu, posz_en) * rozn_obsadzen;
5183 // std::cerr<<"\nrored = "<<rored_posz(E, E0, srednia_masa_el, srednia_masa_dziu, posz_en);
5184 return wynik * struktura::pi / (struktura::c * n_r * struktura::eps0 * E) / struktura::przelm * 1e8;
5185}
5186/*****************************************************************************/
5187 double wzmocnienie::spont_od_pary_poziomow(double E, size_t nr_c, int poz_c, size_t nr_v, int poz_v, double polar)
5188{
5189 // std::cout << "\t\tlumi from pair of levels (0) - " << "nr_c: " << nr_c << ", nr_v: " << nr_v << ", poz_c: " <<
5190 // poz_c << ", poz_v: " << poz_v << "\n";
5191 double wynik;
5192 double cos2tet; // zmiana elementu macierzowego z k_prost
5193 struktura * el = pasma->pasmo_przew[nr_c];
5194 struktura * dziu = pasma->pasmo_wal[nr_v];
5195 /* std::clog<<"E = "<<E<<"\n";
5196 std::clog<<"pasmo c "<<nr_c<<" pasmo v "<<nr_v<<"\n";
5197 std::clog<<"poziomy "<<poz_c<<" i "<<poz_v<<"\n";*/
5198 double Eg; // lokalna przerwa energetyczna
5199 double E0 = Egcv_T[nr_v] - pasma->Egcc[nr_c] + el->rozwiazania[poz_c].poziom +
5200 dziu->rozwiazania[poz_v].poziom; // energia potencjalna + energia prostopadla
5201 // std::cerr<<"spont: poziom c = "<<poz_c<<" poziom v = "<<poz_v<<" E0 = "<<E0<<"\n";
5202 double przekr_w_war;
5203 std::vector<double> calki_kawalki;
5204 double srednia_masa_el = el->rozwiazania[poz_c].prawdopodobienstwa[0] * el->lewa.masa_r;
5205 double srednia_masa_dziu = dziu->rozwiazania[poz_v].prawdopodobienstwa[0] * dziu->lewa.masa_r;
5206 int i;
5207 for(i = 0; i <= (int)el->kawalki.size() - 1; i++)
5208 {
5209 srednia_masa_el += el->rozwiazania[poz_c].prawdopodobienstwa[i + 1] * el->kawalki[i].masa_r;
5210 srednia_masa_dziu += dziu->rozwiazania[poz_v].prawdopodobienstwa[i + 1] * dziu->kawalki[i].masa_r;
5211 }
5212 int ost_ind = el->kawalki.size() + 1;
5213 srednia_masa_el += el->rozwiazania[poz_c].prawdopodobienstwa[ost_ind] * el->prawa.masa_r;
5214 srednia_masa_dziu += dziu->rozwiazania[poz_v].prawdopodobienstwa[ost_ind] * dziu->prawa.masa_r;
5215 double mnoznik_pol;
5216
5217 // double minimalna_przerwa = Egcv_T[nr_v] - pasma->Egcc[nr_c] + el->dol + dziu->dol;
5218 // double min_E0 = Egcv_T[nr_v] - pasma->Egcc[nr_c] + el->rozwiazania[0].poziom + dziu->rozwiazania[0].poziom;
5219 double posz_en; // = posz_z_br(nr_c, poz_c, nr_v, poz_v); // LUKASZ dodalem komentarz
5220 // MD switch-case (lepsze niż if)
5221 switch(wersja)
5222 {
5223 case Z_CHROPOWATOSCIA:
5224 posz_en = posz_z_chrop(nr_c, poz_c, nr_v, poz_v);
5225 break;
5226 case Z_POSZERZENIEM:
5227 posz_en = posz_z_br(nr_c, poz_c, nr_v, poz_v);
5228 break;
5229 }
5230
5231 // double erf_dor = erf_dorored(E, E0, posz_en);
5232 double srednie_k_zeznakiem =
5233 (E - E0 > 0) ? kodE(E - E0, srednia_masa_el, srednia_masa_dziu) : -kodE(E0 - E, srednia_masa_el, srednia_masa_dziu);
5234 Eg = Egcv_T[nr_v] - pasma->Egcc[nr_c];
5235 double sr_E_E0_dod = posz_en / (sqrt(2 * struktura::pi)) * exp(-(E - E0) * (E - E0) / (2 * posz_en * posz_en)) +
5236 (E - E0) * erf_dorored(E, E0, posz_en); // srednia energia kinetyczna w plaszczyznie
5237 // std::clog<<(E-E0)<<" "<<sr_E_E0_dod<<"\n";
5238 // cos2tet= (E0>Eg && E > E0)?(E0-Eg)/(E-Eg):1.0;
5239 // cos2tet = 1.0; // na chwile
5240 cos2tet = (E0 > Eg) ? (E0 - Eg) / (sr_E_E0_dod + E0 - Eg) : 1.;
5241 // std::clog<<"\ncos2tet = "<<cos2tet;
5242 calki_kawalki = (*(pasma->calki_przekrycia_kawalki[nr_c][nr_v]))[poz_c][poz_v];
5243 przekr_w_war = calki_kawalki[0];
5244 mnoznik_pol = (dziu->typ == struktura::hh) ? (1 + cos2tet + polar * (1 - 3 * cos2tet)) / 2
5245 : (5 - 3 * cos2tet + 3 * polar * (3 - cos2tet)) / 6;
5246 // if (dziu->typ == struktura::hh) std::cout << "\t\tlumi from pair of levels (1) - struktura: hh\n"; // LUKASZ
5247 // else std::cout << "\t\tlumi from pair of levels (1) - struktura: lh\n"; // LUKASZ dodalem linie
5248 wynik = sqrt(pasma->el_mac[0] * mnoznik_pol) * przekr_w_war;
5249 for(i = 0; i <= (int)el->kawalki.size() - 1; i++)
5250 {
5251 // std::cerr<<"kawalek "<<i<<"\n";
5252 Eg = Egcv_T[nr_v] - pasma->Egcc[nr_c] + el->kawalki[i].y_pocz +
5253 dziu->kawalki[i].y_pocz; // y_pocz na szybko, moze co innego powinno byc
5254 // cos2tet= (E0>Eg && E > E0)?(E0-Eg)/(E-Eg):1.0;
5255 // cos2tet = 1.0; // na chwile
5256 // cos2tet= (E0>Eg && E > E0)?(E0-Eg)/(sr_E_E0_dod + E0-Eg):1.0;
5257 cos2tet = (E0 > Eg) ? (E0 - Eg) / (sr_E_E0_dod + E0 - Eg) : 1.;
5258 // std::clog<<"\ncos2tet = "<<cos2tet;
5259 mnoznik_pol = (dziu->typ == struktura::hh) ? (1 + cos2tet + polar * (1 - 3 * cos2tet)) / 2
5260 : (5 - 3 * cos2tet + 3 * polar * (3 - cos2tet)) / 6;
5261 // if (dziu->typ == struktura::hh) std::cout << "\t\tlumi from pair of levels (for) - struktura: hh\n"; // LUKASZ
5262 // else std::cout << "\t\tlumi from pair of levels (for) - struktura: lh\n"; // LUKASZ dodalem linie
5263 // std::cerr<<" mnoz_pol = "<<mnoznik_pol<<" cos2tet = "<<cos2tet<<"\n";
5264 przekr_w_war = calki_kawalki[i + 1];
5265 wynik += sqrt(pasma->el_mac[i + 1] * mnoznik_pol) * przekr_w_war;
5266 // std::clog<<"warstwa "<<i<<" przekrycie = "<<przekr_w_war<<"\n";
5267 }
5268 przekr_w_war = calki_kawalki.back();
5269 double energia_elektronu = el->rozwiazania[poz_c].poziom +
5270 srednie_k_zeznakiem * abs(srednie_k_zeznakiem) / (2 * srednia_masa_el); // abs, zeby miec znak, i energie ponizej E0
5271 double energia_dziury =
5272 dziu->rozwiazania[poz_v].poziom + srednie_k_zeznakiem * abs(srednie_k_zeznakiem) / (2 * srednia_masa_dziu);
5273 double obsadzenia = fc(energia_elektronu - pasma->Egcc[nr_c]) * (1 - fv(-energia_dziury + Egcv_T[0] - Egcv_T[nr_v]));
5274 // std::cerr<<"\nE = "<<E<<"\tnr_v = "<<nr_v<<"\ten_el = "<<energia_elektronu<<" en_dziu = "<<energia_dziury<<"
5275 // obsadz el = "<<fc(energia_elektronu - pasma->Egcc[nr_c])<<" obsadz dziu = "<<(1 - fv(-energia_dziury +
5276 // pasma->Egcv[0] - pasma->Egcv[nr_v]))<<" obsadzenia = "<<obsadzenia<<" przesuniecie w fv "<<(pasma->Egcv[0] -
5277 // pasma->Egcv[nr_v])<<"\n";
5278 Eg = Egcv_T[nr_v] - pasma->Egcc[nr_c];
5279 // cos2tet= (E0>Eg && E > E0)?(E0-Eg)/(E-Eg):1.0;
5280 // cos2tet = 1.0; // na chwile
5281 // cos2tet= (E0>Eg && E > E0)?(E0-Eg)/(sr_E_E0_dod + E0-Eg);//:1.0;
5282 cos2tet = (E0 > Eg) ? (E0 - Eg) / (sr_E_E0_dod + E0 - Eg) : 1.;
5283 mnoznik_pol = (dziu->typ == struktura::hh) ? (1 + cos2tet + polar * (1 - 3 * cos2tet)) / 2
5284 : (5 - 3 * cos2tet + 3 * polar * (3 - cos2tet)) / 6;
5285 // if (dziu->typ == struktura::hh) std::cout << "\t\tlumi from pair of levels (2) - struktura: hh\n"; // LUKASZ
5286 // else std::cout << "\t\tlumi from pair of levels (2) - struktura: lh\n"; // LUKASZ dodalem linie
5287 wynik += sqrt(pasma->el_mac.back() * mnoznik_pol) * przekr_w_war;
5288 wynik *= wynik; // dopiero teraz kwadrat modulu
5289 // double posz_en = 2*(E0 - minimalna_przerwa)*pasma->chrop; // oszacowanie rozmycia poziomow z powodu chropowatosci
5290 wynik *= rored_posz(E, E0, srednia_masa_el, srednia_masa_dziu, posz_en) * obsadzenia;
5291 // std::cerr<<"typ_"<<dziu->typ<<" "<<E<<" "<<fc(energia_elektronu - pasma->Egcc[nr_c])<<" "<<(1 - fv(-energia_dziury
5292 // + pasma->Egcv[0] - pasma->Egcv[nr_v]))<<"\n"; std::cerr<<"\nrored = "<<rored_posz(E, E0, srednia_masa_el,
5293 // srednia_masa_dziu, posz_en);
5294 return wynik * E * E * n_r / (struktura::c * struktura::c * struktura::c * struktura::pi * struktura::eps0) /
5295 (struktura::przelm * struktura::przelm * struktura::przelm) * 1e24 / struktura::przels * 1e12; // w 1/(cm^3 s)
5296}
5297/*****************************************************************************/
5298 double wzmocnienie::wzmocnienie_calk_bez_splotu(double E, double polar)
5299{
5300 double wynik = 0.;
5301 for(int nr_c = 0; nr_c <= (int)pasma->pasmo_przew.size() - 1; nr_c++)
5302 for(int nr_v = 0; nr_v <= (int)pasma->pasmo_wal.size() - 1; nr_v++)
5303 wynik += wzmocnienie_od_pary_pasm(E, nr_c, nr_v, polar);
5304 return wynik;
5305}
5306/*****************************************************************************/
5307 double wzmocnienie::wzmocnienie_calk_bez_splotu_L(double lambda, double polar)
5308{
5309 double wynik = 0.;
5310 for(int nr_c = 0; nr_c <= (int) pasma->pasmo_przew.size() - 1; nr_c++)
5311 for(int nr_v = 0; nr_v <= (int) pasma->pasmo_wal.size() - 1; nr_v++)
5312 wynik += wzmocnienie_od_pary_pasm(nm_to_eV(lambda), nr_c, nr_v, polar);
5313 return wynik;
5314}
5315/*****************************************************************************/
5316 void wzmocnienie::profil_wzmocnienia_bez_splotu_dopliku(std::ofstream & plik, double pocz, double kon, double krok)
5317{
5318 double wynikTE, wynikTM;
5319 for(double E = pocz; E <= kon; E += krok)
5320 {
5321 wynikTE = 0; wynikTM = 0;
5322 for(int nr_c = 0; nr_c <= (int)pasma->pasmo_przew.size() - 1; nr_c++)
5323 for(int nr_v = 0; nr_v <= (int)pasma->pasmo_wal.size() - 1; nr_v++)
5324 {
5325 wynikTE += wzmocnienie_od_pary_pasm(E, nr_c, nr_v, 0.);
5326 wynikTE += wzmocnienie_od_pary_pasm(E, nr_c, nr_v, 1.);
5327 std::cerr<<E<<" "<<wynikTE<<" "<<wynikTM<<"\n";
5328 }
5329
5330 plik << E << " " << wynikTE << " " << wynikTM << "\n";
5331 }
5332}
5333/*****************************************************************************/
5334 void wzmocnienie::profil_wzmocnienia_ze_splotem_dopliku(std::ofstream & plik, double pocz, double kon, double krok,
5335 double b)
5336{
5337 for(double E = pocz; E <= kon; E += krok)
5338 {
5339 plik << E << " " << wzmocnienie_calk_ze_splotem(E, b, 0.) << " " << wzmocnienie_calk_ze_splotem(E, b, 1.) << "\n";
5340 }
5341}
5342/*****************************************************************************/
5343 double wzmocnienie::lumin(double E, double polar)
5344{
5345 double wynikTE = 0., wynikTM = 0.;
5346 for (int nr_c = 0; nr_c <= (int)pasma->pasmo_przew.size() - 1; nr_c++)
5347 for (int nr_v = 0; nr_v <= (int)pasma->pasmo_wal.size() - 1; nr_v++)
5348 {
5349 wynikTE += spont_od_pary_pasm(E, nr_c, nr_v, 0.0);
5350 wynikTM += spont_od_pary_pasm(E, nr_c, nr_v, 1.0);
5351 }
5352 if (polar == 0.)
5353 return wynikTE;
5354 else if (polar == 1.)
5355 return wynikTM;
5356 else
5357 return (wynikTE + wynikTM);
5358}
5359/*****************************************************************************/
5360 void wzmocnienie::profil_lumin_dopliku(std::ofstream & plik, double pocz, double kon, double krok)
5361{
5362 // double wynik;
5363 /*
5364 std::vector<double> wklady;
5365 wklady.resize(pasmo_wal.size());
5366 */
5367 double wynikTE, wynikTM;
5368 for(double E = pocz; E <= kon; E += krok)
5369 {
5370 plik << E;
5371 wynikTE = 0.0;
5372 wynikTM = 0.0;
5373 for(int nr_c = 0; nr_c <= (int)pasma->pasmo_przew.size() - 1; nr_c++)
5374 for(int nr_v = 0; nr_v <= (int)pasma->pasmo_wal.size() - 1; nr_v++)
5375 {
5376 wynikTE += spont_od_pary_pasm(E, nr_c, nr_v, 0.0);
5377 wynikTM += spont_od_pary_pasm(E, nr_c, nr_v, 1.0);
5378 }
5379 plik << "\t" << wynikTE << " " << wynikTM;
5380 // plik<<E<<" "<<wynik<<"\n";
5381 plik << std::endl;
5382 }
5383}
5384/*****************************************************************************/
5385 void wzmocnienie::profil_lumin_dopliku_L(std::ofstream & plik, double pocz, double kon, double krok) // dodalem metode LUKASZ LUKI 5.09.2023
5386{
5387 double wynikTE, wynikTM;
5388 for(double L = pocz; L <= kon; L += krok)
5389 {
5390 wynikTE = 0.0;
5391 wynikTM = 0.0;
5392 for(int nr_c = 0; nr_c <= (int) pasma->pasmo_przew.size() - 1; nr_c++)
5393 for(int nr_v = 0; nr_v <= (int) pasma->pasmo_wal.size() - 1; nr_v++)
5394 {
5395 wynikTE += spont_od_pary_pasm(nm_to_eV(L), nr_c, nr_v, 0.0);
5396 wynikTM += spont_od_pary_pasm(nm_to_eV(L), nr_c, nr_v, 1.0);
5397 }
5398 plik<<L<<" "<<wynikTE<<" "<<wynikTM<<"\n";
5399 }
5400}
5401/*****************************************************************************/
5402 double wzmocnienie::moc_lumin()
5403{
5404 struktura * el = pasma->pasmo_przew[0];
5405 struktura * dziu = pasma->pasmo_wal[0];
5406 double min_E0 = Egcv_T[0] - pasma->Egcc[0] + el->rozwiazania[0].poziom + dziu->rozwiazania[0].poziom;
5407 double lok_min_E0;
5408 for(int nr_c = 0; nr_c <= (int)pasma->pasmo_przew.size() - 1; nr_c++)
5409 for(int nr_v = 0; nr_v <= (int)pasma->pasmo_wal.size() - 1; nr_v++)
5410 {
5411 lok_min_E0 = Egcv_T[nr_v] - pasma->Egcc[nr_c] + el->rozwiazania[0].poziom + dziu->rozwiazania[0].poziom;
5412 min_E0 = (min_E0 < lok_min_E0) ? min_E0 : lok_min_E0;
5413 }
5414 double minimalna_przerwa = pasma->min_przerwa_energetyczna();
5415 double posz_en = 2 * (min_E0 - minimalna_przerwa) * pasma->chrop;
5416 double pocz = min_E0 - 2 * posz_en;
5417 double kon = min_E0 + 6 * struktura::kB * T;
5418 kon = (pocz < kon) ? kon : pocz + 2 * struktura::kB * T;
5419#ifdef LOGUJ // MD
5420 std::clog << "\nW mocy. pocz = " << pocz << " kon = " << kon << "\n";
5421#endif
5422 double krok = struktura::kB * T / 30;
5423 double wynik = 0;
5424 for(double E = pocz; E <= kon; E += krok)
5425 {
5426 for(int nr_c = 0; nr_c <= (int)pasma->pasmo_przew.size() - 1; nr_c++)
5427 for(int nr_v = 0; nr_v <= (int)pasma->pasmo_wal.size() - 1; nr_v++)
5428 {
5429 wynik += spont_od_pary_pasm(E, nr_c, nr_v, 0.0); // 0.0 dla TE
5430 }
5431 }
5432 return wynik * krok;
5433}
5434/*****************************************************************************/
5435 double wzmocnienie::fc(double E)
5436{
5437 double arg = (E - qFlc) / (struktura::kB * T);
5438 return 1 / (1 + exp(arg));
5439}
5440/*****************************************************************************/
5441 double wzmocnienie::fv(double E)
5442{
5443 double arg = (E - qFlv) / (struktura::kB * T);
5444 return 1 / (1 + exp(arg));
5445}
5446/*****************************************************************************/
5447 double wzmocnienie::przel_gest_z_cm2(double gest_w_cm2) // gestosc powierzchniowa
5448{
5449 return gest_w_cm2 * 1e-16 * struktura::przelm * struktura::przelm;
5450}
5451/*****************************************************************************/
5452 double wzmocnienie::przel_gest_na_cm2(double gest_w_wew) // gestosc powierzchniowa
5453{
5454 return gest_w_wew / (1e-16 * struktura::przelm * struktura::przelm);
5455}
5456/*****************************************************************************/
5457 double wzmocnienie::posz_z_chrop(size_t nr_c, int poz_c, size_t nr_v, int poz_v)
5458{
5459 double srdno;
5460 double E0;
5461 double Eklok;
5462 double srednia_ekz_el = 0;
5463 double srednia_ekz_dziu = 0;
5464 double grub;
5465 struktura * el = pasma->pasmo_przew[nr_c];
5466 struktura * dziu = pasma->pasmo_wal[nr_v];
5467 double posz_en = 0;
5468 for(int i = 0; i <= (int)el->kawalki.size() - 1; i++)
5469 {
5470 // std::cerr<<"\n prawd "<<i<<" = "<<el->rozwiazania[poz_c].prawdopodobienstwa[i+1];
5471 srdno = (el->kawalki[i].y_pocz + el->kawalki[i].y_kon) / 2;
5472 grub = (el->kawalki[i].x_kon - el->kawalki[i].x_pocz);
5473 E0 = el->rozwiazania[poz_c].poziom;
5474 Eklok = E0 - srdno;
5475 srednia_ekz_el =
5476 (Eklok > 0) ? el->rozwiazania[poz_c].prawdopodobienstwa[i + 1] * Eklok : 0; // srednia energia kinetyczna po z
5477 srdno = (dziu->kawalki[i].y_pocz + dziu->kawalki[i].y_kon) / 2;
5478 E0 = dziu->rozwiazania[poz_v].poziom;
5479 Eklok = E0 - srdno;
5480 srednia_ekz_dziu = (Eklok > 0) ? dziu->rozwiazania[poz_v].prawdopodobienstwa[i + 1] * Eklok : 0;
5481 posz_en += 2 * (srednia_ekz_el + srednia_ekz_dziu) * pasma->chrop / grub;
5482 // std::clog<<"\npoz_c = "<<poz_c<<" poz_v = "<<poz_v<<"\tn_srdno = "<<srdno<<"\tE0 = "<<E0<<"\tsrEk =
5483 // "<<srednia_ekz_el<<"\tposze_en = "<<posz_en; srednia_ekz_dziu +=
5484 // dziu->rozwiazania[poz_v].prawdopodobienstwa[i+1]*dziu->kawalki[i].masa_r;
5485 }
5486 return posz_en;
5487}
5488/*****************************************************************************/
5489 double wzmocnienie::posz_z_br(size_t nr_c, int poz_c, size_t nr_v, int poz_v)
5490{
5491 // std::clog<<"Jestem w posz_z_br\n";
5492 struktura * el1 = pasma->pasmo_przew[nr_c];
5493 struktura * el2 = pasma->pasmo_przew_mod[nr_c];
5494 struktura * dziu1 = pasma->pasmo_wal[nr_v];
5495 struktura * dziu2 = pasma->pasmo_wal_mod[nr_v];
5496 if(el2 == NULL || dziu2 == NULL)
5497 {
5498 Error err; // MD
5499 err << "\nNie ma drugiej struktury!\n";
5500 throw err;
5501 }
5502 double dEe;
5503 double dEd;
5504 int rozmmodel = el2->rozwiazania.size();
5505 int rozmmoddziu = dziu2->rozwiazania.size();
5506 if(poz_c > rozmmodel - 1) // brakuje poziomu mod
5507 {
5508 dEe = el1->rozwiazania[rozmmodel - 1].poziom - el2->rozwiazania[rozmmodel - 1].poziom;
5509 }
5510 else
5511 dEe = el1->rozwiazania[poz_c].poziom - el2->rozwiazania[poz_c].poziom;
5512 if(poz_v > rozmmoddziu - 1) // brakuje poziomu mod
5513 {
5514 dEd = dziu1->rozwiazania[rozmmoddziu - 1].poziom - dziu2->rozwiazania[rozmmoddziu - 1].poziom;
5515 // dEd = 0.0;
5516 }
5517 else
5518 dEd = dziu1->rozwiazania[poz_v].poziom - dziu2->rozwiazania[poz_v].poziom;
5519 double dE = dEe + dEd;
5520 dE = (dE >= 0) ? dE : -dE;
5521 // std::clog<<"pozc,v ="<<poz_c<<", "<<poz_v<<"\tdE = "<<dE<<"\n";
5522 return dE * pasma->broad; // tutaj chrop oznacza mnoznik dla roznicy miedzy pasmami a pasmami_mod
5523}
5524
5525} // namespace kubly