📄 d9r2f.java
字号:
public class d9r2F
{
double lift_chisq;
void lfit(double x[], double y[], double sig[], int ndata, double a[], int ma,
int lista[], int mfit, double covar[], int ncvm, double chisq)
{
int i, j, k, kk, ihit;
double ym, sig2i ,wt, sum;
double beta[] = new double[51];
double afunc[] = new double[51];
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 set in lista");
System.exit(1);
}
}
}
if (kk != (ma + 1))
{
System.out.println(" improper set in lista");
System.exit(1);
}
for (j = 1; j <= mfit; j++)
{
for (k = 1; k <= mfit; k++)
{
covar[(j - 1) * ma + k] = 0.0;
}
beta[j] = 0.0;
}
for (i = 1; i <= ndata; i++)
{
funcs(x[i], afunc, ma);
ym = y[i];
if (mfit < ma)
{
for (j = mfit + 1; j <= ma; j++)
{
ym = ym - a[lista[j]] * afunc[lista[j]];
}
}
sig2i = 1.0 / (sig[i] * sig[i]);
for (j = 1; j <= mfit; j++)
{
wt = afunc[lista[j]] * sig2i;
for (k = 1; k <= j; k++)
{
covar[(j - 1) * ma + k] = covar[(j - 1) * ma + k] + wt * afunc[lista[k]];
}
beta[j] = beta[j] + ym * wt;
}
}
if (mfit > 1)
{
for (j = 2; j <= mfit; j++)
{
for (k = 1; k <= j - 1; k++)
{
covar[(k - 1) * ma + j] = covar[(j - 1) * ma + k];
}
}
}
gaussj(covar, mfit, ma, beta);
for (j = 1; j <= mfit; j++)
{
a[lista[j]] = beta[j];
}
chisq = 0.0;
for (i = 1; i <= ndata; i++)
{
funcs(x[i], afunc, ma);
sum = 0.0;
for (j = 1; j <= ma; j++)
{
sum = sum + a[j] * afunc[j];
}
chisq = chisq + ((y[i] - sum) / sig[i]) * ((y[i] - sum) / sig[i]);
}
covsrt(covar, ncvm, ma, lista, mfit);
lift_chisq = chisq;
}
void funcs(double x, double afunc[], int ma)
{
afunc[1] = 1.0;
for(int i = 2; i <= ma; i++)
{
afunc[i] = x * afunc[i - 1];
}
}
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;
}
}
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;
}
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;
}
}
}
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -