--------------------------------------------------------------------------------------------------------------------
Expand|Select|Wrap|Line Numbers
- double a, b, c, d;
- int n_max, n; //n_max = number of iterations & n = running index
- cout << endl << "enter guess for a ... " << endl;
- cin >> a;
- cout << endl << "enter guess for b ... " << endl;
- cin >> b;
- cout << endl << "enter guess for c ... " << endl;
- cin >> c;
- cout << endl << "enter maximum number of iterations ... " << endl;
- cin >> n_max;
- long double K_old[3][1], K_new[3][1], K_min[3][1];
- K_old [1][1] = a;
- K_old [2][1] = b;
- K_old [3][1] = c;
- double chi_min;
- chi_min = 2000;
- long double chi[n_max];
- for(n = 0; n <= n_max; n++)
- {
- long double Chi, Chi_prime[1][3], Chi_prime_prime[3][3][n_max];
- double X, i[9096];
- int x;
- Chi = 0;
- Chi_prime_prime[1][1][n] = 0;
- Chi_prime_prime[1][2][n] = 0;
- Chi_prime_prime[1][3][n] = 0;
- Chi_prime_prime[2][1][n] = 0;
- Chi_prime_prime[2][2][n] = 0;
- Chi_prime_prime[2][3][n] = 0;
- Chi_prime_prime[3][1][n] = 0;
- Chi_prime_prime[3][2][n] = 0;
- Chi_prime_prime[3][3][n] = 0;
- X = 2 * S * S * -1;
- for (j = 0; j < m; j++) //calculating grad's and concavity
- {
- i[j] = Imax * exp((lambda[j] - average)/X) - (a * a * lambda[j] + b * lambda[j] + c);
- Chi = Chi + pow(((I[j] - i[j])/S), 2);
- Chi_prime[1][1] = Chi_prime[1][1] + 2*((I[j] - i[j])/S)*((2 * a * lambda[j])/S); //derivative w.r.t. a
- Chi_prime[2][1] = Chi_prime[2][1] + 2*((I[j] - i[j])/S)*(lambda[j]/S); //derivative b
- Chi_prime[3][1] = Chi_prime[3][1] + 2*((I[j] - i[j])/S)*(1/S); //derivative c
- Chi_prime_prime[1][1][n] = Chi_prime_prime[1][1][n] + ((I[j] - i[j])/S)*((4 * lambda[j])/S) + ((4 * a * lambda[j])/S)*((2 * a * lambda[j])/S);
- Chi_prime_prime[1][2][n] = Chi_prime_prime[1][2][n] + (((4 * a * lambda[j])/S)*(lambda[j]/S));
- Chi_prime_prime[1][3][n] = Chi_prime_prime[1][3][n] + ((4 * a * lambda[j])/S)*(1/S);
- Chi_prime_prime[2][1][n] = Chi_prime_prime[2][1][n] + ((2 * lambda[j])/S)*((2 * a * lambda[j])/S);
- Chi_prime_prime[2][2][n] = Chi_prime_prime[2][2][n] + (((2 * lambda[j])/S)*(lambda[j]/S));
- Chi_prime_prime[2][3][n] = Chi_prime_prime[2][3][n] + (((2 * lambda[j])/S)*(1/S));
- Chi_prime_prime[3][1][n] = Chi_prime_prime[3][1][n] + (2/S)*((2 * a * lambda[j])/S);
- Chi_prime_prime[3][2][n] = Chi_prime_prime[3][2][n] + ((2/S)*(lambda[j]/S));
- Chi_prime_prime[3][3][n] = Chi_prime_prime[3][3][n] + ((2/S)*(1/S));
- }
- chi[n] = Chi;
- cout << "chi squared = " << chi[n] << " ";
- prt4 << "chi squared = " << chi[n] << " ";
- jacobian << "n = " << n << endl;
- jacobian << Chi_prime_prime[1][1][n] << " ";
- jacobian << Chi_prime_prime[1][2][n] << " ";
- jacobian << Chi_prime_prime[1][3][n] << endl;
- jacobian << Chi_prime_prime[2][1][n] << " ";
- jacobian << Chi_prime_prime[2][2][n] << " ";
- jacobian << Chi_prime_prime[2][3][n] << endl;
- jacobian << Chi_prime_prime[3][1][n] << " ";
- jacobian << Chi_prime_prime[3][2][n] << " ";
- jacobian << Chi_prime_prime[3][3][n] << endl << endl;
- long double det_J = Chi_prime_prime[1][1][n]*(Chi_prime_prime[2][2][n]*Chi_prime_prime[3][3][n]
- - (Chi_prime_prime[2][3][n]*Chi_prime_prime[3][2][n])) - Chi_prime_prime[1][2][n]
- *(Chi_prime_prime[2][1][n]*Chi_prime_prime[3][3][n] - (Chi_prime_prime[2][3][n]*
- Chi_prime_prime[3][1][n])) + Chi_prime_prime[1][3][n]*(Chi_prime_prime[2][1][n]*
- Chi_prime_prime[3][2][n] - (Chi_prime_prime[2][2][n]*Chi_prime_prime[3][1][n]));
- long double inv_J[3][3][n];
- inv_J[1][1][n] = (1/det_J)*Chi_prime_prime[1][1][n];
- inv_J[1][2][n] = (1/det_J)*Chi_prime_prime[1][2][n];
- inv_J[1][3][n] = (1/det_J)*Chi_prime_prime[1][3][n];
- inv_J[2][1][n] = (1/det_J)*Chi_prime_prime[2][1][n];
- inv_J[2][2][n] = (1/det_J)*Chi_prime_prime[2][2][n];
- inv_J[2][3][n] = (1/det_J)*Chi_prime_prime[2][3][n];
- inv_J[3][1][n] = (1/det_J)*Chi_prime_prime[3][1][n];
- inv_J[3][2][n] = (1/det_J)*Chi_prime_prime[3][2][n];
- inv_J[3][3][n] = (1/det_J)*Chi_prime_prime[3][3][n];
- inverse << "n = " << n << endl;
- inverse << inv_J[1][1][n] << " ";
- inverse << inv_J[1][2][n] << " ";
- inverse << inv_J[1][3][n] << endl;
- inverse << inv_J[2][1][n] << " ";
- inverse << inv_J[2][2][n] << " ";
- inverse << inv_J[2][3][n] << endl;
- inverse << inv_J[3][1][n] << " ";
- inverse << inv_J[3][2][n] << " ";
- inverse << inv_J[3][3][n] << endl << endl;
- long double D[3][1];
- for (0 < x; x <= 3; x++)
- D[x][1] = 0;
- for (0 < x; x <= 3; x++)
- D[1][1] = D[1][1] + inv_J[1][x][n]*Chi_prime[1][1];
- D[2][1] = D[2][1] + inv_J[2][x][n]*Chi_prime[2][1];
- D[3][1] = D[3][1] + inv_J[3][x][n]*Chi_prime[3][1];
- prt4 << "D[1][1] = " << D[1][1] << ", D[2][1] = " << D[2][1];
- prt4 << ", D[3][1] = " << D[3][1] << endl << endl;
- cout << "D[1][1] = " << D[1][1] << ", D[2][1] = " << D[2][1];
- cout << ", D[3][1] = " << D[3][1] << endl;
- K_new[1][1] = K_old[1][1] - D[1][1];
- K_new[2][1] = K_old[2][1] - D[2][1];
- K_new[3][1] = K_old[3][1] - D[3][1];
- prt4 << "K[1][1] = " << K_new[1][1] << ", K[2][1] = " << K_new[2][1];
- prt4 << ", K[3][1] = " << K_new[3][1] << endl << endl;
- cout << "K[1][1] = " << K_new[1][1] << ", K[2][1] = " << K_new[2][1];
- cout << ", K[3][1] = " << K_new[3][1] << endl;
- if (n > 0)
- if (Chi_prime_prime[1][1][n] > 0 && Chi_prime_prime[1][1][n - 1] > 0)//concavity stays const.
- if (Chi_prime_prime[2][2][n] > 0 && Chi_prime_prime[2][2][n - 1] > 0)
- if (Chi_prime_prime[3][3][n] > 0 && Chi_prime_prime[3][3][n - 1] > 0)
- {
- if (chi_min - chi [n] > 0)
- {
- chi_min = chi[n];
- K_min[1][1] = K_new[1][1];
- K_min[2][1] = K_new[2][1];
- K_min[3][1] = K_new[3][1];
- }
- else {
- chi_min = chi_min;
- K_min[1][1] = K_min[1][1];
- K_min[2][1] = K_min[2][1];
- K_min[3][1] = K_min[3][1];
- }
- }
- else {
- prt4 << "point of inflection passed or in vicinity of local max" << endl << endl;
- cout << "minimised chi = " << chi_min << " ";
- cout << "corresponding K = " << K_min[1][1];
- cout << ", " << K_min[2][1];
- cout << ", " << K_min[3][1] << endl << endl;
- prt << "minimised chi = " << chi_min;
- prt << " corresponding K = " << K_min[1][1];
- prt << ", " << K_min[2][1];
- prt << ", " << K_min[3][1] << endl << endl;
- for (j = 0; j < m; j++)
- {
- //cout << "i[" << j + 1 << "] = " << i[j] << endl;
- //intensity << "i[" << j + 1 << "] = " << i[j] << endl;
- intensity << i[j] << endl;
- }
- system ("pause");
- exit(1);
- }
- if (n > 2) //mimimisation criteria
- if ((chi[n] = chi_min) && (fabs(chi[n - 1] - chi_min)) < 0.1 && (fabs(chi[n - 2] - chi_min)) < 0.1)
- {
- cout << "minimised chi = " << chi_min << " " << "corresponding K = ";
- cout << K_min[1][1] << ", " << K_min[2][1] << ", " << K_min[3][1] << endl << endl;
- prt << "minimised chi = " << chi_min << " " << "corresponding K = ";
- prt << K_min[1][1] << ", " << K_min[2][1] << ", " << K_min[3][1] << endl << endl;
- for (j = 0; j < m; j++)
- {
- //cout << "i[" << j + 1 << "] = " << i[j] << endl;
- //intensity << "i[" << j + 1 << "] = " << i[j] << endl;
- intensity << i[j] << endl;
- }
- system ("pause");
- exit(1);
- }
- K_old[1][1] = K_new[1][1];
- K_old[2][1] = K_new[2][1];
- K_old[3][1] = K_new[3][1];
- }
- cout << "minimised chi = " << chi_min << " " << "corresponding K = ";
- cout << K_min[1][1] << ", " << K_min[2][1] << ", " << K_min[3][1] << endl;
- prt << "minimised chi = " << chi_min << " " << "corresponding K = ";
- prt << K_min[1][1] << ", " << K_min[2][1] << ", " << K_min[3][1] << endl;
- prt.close();
- prt4.close();
- system ("pause");
- }