📄 d9r8af.java
字号:
public class d9r8aF
{
double fgauss_y;
void fgauss(double x, double a[], double y, double dyda[], int na)
{
y = 0.0;
for (int i = 1; i <= na - 1; i = i + 3)
{
double arg = (x - a[i + 1]) / a[i + 2];
double ex = Math.exp(-(arg * arg));
double fac = a[i] * ex * 2.0 * arg;
y = y + a[i] * ex;
dyda[i] = ex;
dyda[i + 1] = fac / a[i + 2];
dyda[i + 2] = fac * arg / a[i + 2];
}
fgauss_y = y;
}
double mrqcof_chisq;
void mrqcof(double x[], double y[], double sig[], int ndata, double a[],
int ma, int lista[], int mfit, double alpha[], double beta[],
int nalp, double chisq)
{
int i, j, k;
double wt, ymod, sig2i, dy;
double dyda[] = new double[21];
ymod = 0;
for (j = 1; j <= mfit; j++)
{
for (k = 1; k <= j; k++)
{
alpha[(j - 1) * nalp + k] = 0.0;
}
beta[j] = 0.0;
}
chisq = 0.0;
for (i = 1; i <= ndata; i++)
{
fgauss(x[i], a, ymod, dyda, ma);
ymod = fgauss_y;
sig2i = 1.0 / (sig[i] * sig[i]);
dy = y[i] - ymod;
for (j = 1; j <= mfit; j++)
{
wt = dyda[lista[j]] * sig2i;
for (k = 1; k <= j; k++)
{
alpha[(j - 1) * nalp + k] = alpha[(j - 1) * nalp + k] + wt * dyda[lista[k]];
}
beta[j] = beta[j] + dy * wt;
}
chisq = chisq + dy * dy * sig2i;
}
for (j = 2; j <= mfit; j++)
{
for (k = 1; k <= j - 1; k++)
{
alpha[(k - 1) * nalp + j] = alpha[(j - 1) * nalp + k];
}
}
mrqcof_chisq = chisq;
}
static int kk, i, j, k, ihit;
static double ochisq;
double mrqmin_chisq;
double mrqmin_alamda;
static double atry[] = new double[21];
static double beta[] = new double[21];
static double da[] = new double[21];
void mrqmin(double x[], double y[], double sig[], int ndata, double a[],
int ma, int lista[], int mfit, double covar[], double alpha[],
int nca, double chisq, double alamda)
{
//double ochisq=0;
if (alamda < 0.0)
{
kk = mfit + 1;
for (j = 1; j <= ma; j++)
{
ihit = 0;
for (k = 1; k <= mfit; k++)
{
if (lista[k] == j)
{
ihit = ihit + 1;
}
}
if (ihit == 0)
{
lista[kk] = j;
kk = kk + 1;
}
else
{
if (ihit > 1)
{
System.out.println ("improper permutation in lista");
System.exit(1);
}
}
}
if (kk != ma + 1)
{
System.out.println ("improper permutation in lista");
System.exit(1);
}
alamda = 0.001;
mrqcof(x, y, sig, ndata, a, ma, lista, mfit, alpha, beta, nca, chisq);
chisq = mrqcof_chisq;
ochisq = chisq;
for (j = 1; j <= ma; j++)
{
atry[j] = a[j];
}
}
for (j = 1; j <= mfit; j++)
{
for (k = 1; k <= mfit; k++)
{
covar[(j - 1) * nca + k] = alpha[(j - 1) * nca + k];
}
covar[(j - 1) * nca + j] = alpha[(j - 1) * nca + j] * (1.0 + alamda);
da[j] = beta[j];
}
gaussj(covar, mfit, nca, da);
if (alamda == 0.0)
{
covsrt(covar, nca, ma, lista, mfit);
for (i = 1; i <= 20; i++)
{
da[i] = 0;
atry[i] = 0;
}
mrqmin_chisq = chisq; mrqmin_alamda = alamda;
return;
}
for (j = 1; j <= mfit; j++)
{
atry[lista[j]] = a[lista[j]] + da[j];
}
mrqcof(x, y, sig, ndata, atry, ma, lista, mfit, covar, da, nca, chisq);
chisq = mrqcof_chisq;
if (chisq < ochisq)
{
alamda = 0.1 * alamda;
ochisq = chisq;
for (j = 1; j <= mfit; j++)
{
for (k = 1; k <= mfit; k++)
{
alpha[(j - 1) * nca + k] = covar[(j - 1) * nca + k];
}
beta[j] = da[j];
int temp = lista[j];
a[temp] = atry[temp];
//a[lista[j]] = atry[lista[j]];
}
}
else
{
alamda = 10.0 * alamda;
chisq = ochisq;
}
for (i=1; i <= 20; i++)
{
da[i] = 0;
atry[i] = 0;
}
mrqmin_chisq = chisq; mrqmin_alamda = alamda;
}
static int iset;
static double gset;
int gasdev_idum;
double gasdev(int idum)
{
double t, v1, v2, fac, r;
if (iset == 0.0)
{
do
{
v1 = 2.0 * ran1(idum) - 1.0;
idum = ran1_idum;
v2 = 2.0 * ran1(idum) - 1.0;
idum = ran1_idum;
r = v1 * v1 + v2 * v2;
}while ((r >= 1.0) || (r == 0.0));
fac = Math.sqrt(-2.0 * Math.log(r) / r);
gset = v1 * fac;
t = v2 * fac;
iset = 1;
gasdev_idum = idum;
return t;
}
else
{
t = gset;
gasdev_idum = idum;
iset = 0;
return t;
}
}
void covsrt(double covar[], int ncvm, int ma, int lista[], int mfit)
{
int i, j;
double swap;
for (j = 1; j <= ma - 1; j++)
{
for (i = j + 1; i <= ma; i++)
{
covar[(i - 1) * ma + j] = 0.0;
}
}
for (i = 1; i <= mfit - 1; i++)
{
for (j = i + 1; j <= mfit; j++)
{
if (lista[j] > lista[i])
{
covar[(lista[j] - 1) * ma + lista[i]] = covar[(i - 1) * ma + j];
}
else
{
covar[(lista[i] - 1) * ma + lista[j]] = covar[(i - 1) * ma + j];
}
}
}
swap = covar[1];;
for (j = 1; j <= ma; j++)
{
covar[j] = covar[(j - 1) * ma + j];
covar[(j - 1) * ma + j] = 0.0;
}
covar[(lista[1] - 1) * ma + lista[1]] = swap;
for (j = 2; j <= mfit; j++)
{
covar[(lista[j] - 1) * ma + lista[j]] = covar[j];
}
for (j = 2; j <= ma; j++)
{
for (i = 1; i <= j - 1; i++)
{
covar[(i - 1) * ma + j] = covar[(j - 1) * ma + i];
}
}
}
void gaussj(double A[], int N, int MA,double B[])
{
int i, j, k, l ,ll, IROW, ICOL;
double BIG,PIVINV,DUM;
int IPIV[] = new int[51];
int INDXR[] = new int[51];
int INDXC[] = new int[51];
ICOL = 0;
IROW = 0;
for (j = 1; j <= N; j++)
{
IPIV[j] = 0;
}
for (i = 1; i <= N; i++)
{
BIG = 0.0;
for (j = 1; j <= N; j++)
{
if(IPIV[j] != 1)
{
for(k = 1; k <= N; k++)
{
if(IPIV[k] == 0)
{
if(Math.abs(A[(j - 1) * MA + k]) >= BIG)
{
BIG = Math.abs(A[(j - 1) * MA + k]);
IROW = j;
ICOL = k;
}
else if(IPIV[k] > 1)
{
System.out.println("Singular matrix");
System.exit(1);
}
}
}
}
}
IPIV[ICOL]=IPIV[ICOL] + 1;
if(IROW != ICOL)
{
for (l = 1; l <= N; l++)
{
DUM = (A[(IROW - 1) * MA + l]);
A[(IROW - 1) * MA + l] = A[(ICOL - 1) * MA + l];
A[(ICOL - 1) * MA + l] = DUM;
}
DUM = B[IROW];
B[IROW] = B[ICOL];
B[ICOL] = DUM;
}
INDXR[i] = IROW;
INDXC[i] = ICOL;
if(A[(ICOL - 1) * MA + ICOL] == 0.0)
{
System.out.println("Singular matrix");
System.exit(1);
}
PIVINV = 1.0 / (A[(ICOL - 1) * MA + ICOL]);
A[(ICOL - 1) * MA + ICOL] = 1.0;
for (l = 1; l <= N; l++)
{
A[(ICOL - 1) * MA + l] = A[(ICOL - 1) * MA + l] * PIVINV;
}
B[ICOL] = B[ICOL] * PIVINV;
for (ll = 1; ll <= N; ll++)
{
if(ll != ICOL)
{
DUM = A[(ll - 1) * MA + ICOL];
A[(ll - 1) * MA + ICOL] = 0.0;
for(l = 1; l <= N; l++)
{
A[(ll - 1) * MA + l] = A[(ll - 1) * MA + l] - A[(ICOL - 1) * MA + l] * DUM;
}
B[ll] = B[ll] - B[ICOL] * DUM;
}
}
}
for (l = N; l >= 1; l--)
{
if(INDXR[l] != INDXC[l])
{
for (k = 1; k <= N; k++)
{
DUM = A[(k - 1) * MA + INDXR[l]];
A[(k - 1) * MA + INDXR[l]] = A[(k - 1) * MA + INDXC[l]];
A[(k - 1) * MA + INDXR[l]] = DUM;
}
}
}
}
int ran1_idum;
static double r[] = new double[98];
static int ix1, ix2, ix3;
static int iff;
double ran1(int idum)
{
double rm1, rm2, t;
int m1, m2, m3, ia1, ia2, ia3, ic1, ic2, ic3, j;
m1 = 259200; ia1 = 7141; ic1 = 54773; rm1 = 0.0000038580247;
m2 = 134456; ia2 = 8121; ic2 = 28411; rm2 = 0.0000074373773;
m3 = 243000; ia3 = 4561; ic3 = 51349;
if ((idum < 0) || (iff == 0))
{
iff = 1;
ix1 = (ic1 - idum) % m1;
ix1 = (ia1 * ix1 + ic1) % m1;
ix2 = ix1 % m2;
ix1 = (ia1 * ix1 + ic1) % m1;
ix3 = ix1 % m3;
for (j = 1; j <= 97; j++)
{
ix1 = (ia1 * ix1 + ic1) % m1;
ix2 = (ia2 * ix2 + ic2) % m2;
r[j] = ((double)ix1 + (double)ix2 * rm2) * rm1;
}
idum = 1;
}
ix1 = (ia1 * ix1 + ic1) % m1;
ix2 = (ia2 * ix2 + ic2) % m2;
ix3 = (ia3 * ix3 + ic3) % m3;
j = 1 + (int)((97 * ix3) / m3);
if ((j > 97) || (j < 1))
{
System.out.println( "abnormal exit");
System.exit(1);
}
t = r[j];
r[j] = ((double)ix1 + (double)ix2 * rm2) * rm1;
ran1_idum = idum;
return t;
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -