Разложение импульса.Разложение ЭДС.Использование регуляризации, страница 14

{

if( t1.y*t3.y >= 0 )s += (t1.y + t3.y)*(t1.x - t3.x)/2.0;else

{

x = -(t3.x - t1.x)/(t3.x - t1.y)*t1.y + t1.x;

s += (t1.y*(t1.x - x) + t3.y*(x - t3.x))/2.0;

}

s /= (cur.x - t3.x);

Coeff[m*2] = s - sum;

Coeff[(m++)*2+1] = cur.x;

sum = s;

s = 0;

epsi *= cthick;

t1 = cur = t3;

}else cur.y = t3.y;

}else

{

if( ProgressFunc != NULL )

if( !(*ProgressFunc)( paContinue, (in-i)*100/in ) )return -10;

if( t1.y*t2.y >= 0 )s += (t1.y + t2.y)*(t1.x - t2.x)/2.0;else

{

x = -(t2.x - t1.x)/(t2.y - t1.y)*t1.y + t1.x;

s += (t1.y*(t1.x - x) + t2.y*(x - t2.x))/2.0;

}

if( i < 0 )break;

t1 = t2;

t2.x = IT[i]; t2.y = II[i--];

}

}

if( fabs(cur.x-t2.x) > epst )

{

s /= (cur.x - t2.x);

Coeff[m*2] = s - sum;

Coeff[(m++)*2+1] = cur.x;

sum = s;

}

Coeff[m*2] = -sum;

Coeff[(m++)*2+1] = t2.x;

return errOk;

}

int CEMFTransform::DelImpulsTime()

{

int res = DecomposeImpuls();

if( res != errOk )return res;

if( !ET || !EMF )return errNoEMFCurve;

while( ET[0] < Coeff[1] && n > 0 )

{

ET++;

EMF++;

n--;

}

return errOk;

}

int CEMFTransform::MakeMatrix( double *M, int *IM, double *HT, int &hn )

{

int i, j, k;

double t,srct;

hn = n+1; // One additional point required

// Make heaviside times

HT[hn-1] = ET[hn-2] - Coeff[2*m-1];

for( i = 0; i < hn-1; i++ )

HT[i] = ET[i] - Coeff[1];

if( ProgressFunc != NULL )

(*ProgressFunc)( paMakeMatrix, 0 );

memset( M, 0, hn*(hn+1)/2*sizeof(double) );

IM[0] = 0;

for( j = 0; j < hn-1; j++ )

{

if( ProgressFunc != NULL )

if( !(*ProgressFunc)( paContinue, j*100/(hn-2) ) )return errUserBreak;

srct = HT[j] + Coeff[1];

M[IM[j]] = Coeff[0];

k = j+1;

for( i = 1; i < m; i++ )

{

t = srct - Coeff[2*i+1];

while( HT[k] < t )k++;

k--;

t = (t - HT[k])/(HT[k+1] - HT[k]);

M[IM[j]+k-j] += Coeff[2*i]*(1-t);

M[IM[j]+k-j+1] += Coeff[2*i]*t;

}

IM[j+1] = IM[j] + k - j + 2;

}

IM[hn] = IM[hn-1] + 1;

// Make last equation

M[IM[hn-1]] = M[IM[hn-2]] + pow(HT[hn-2],2.5)*pow(HT[hn-1],-2.5)*M[IM[hn-2]+1];

EMF[hn-1] = EMF[hn-2]*pow(HT[hn-2],2.5)*pow(HT[hn-1],-2.5);

// -------------------------return errOk;

}

void CEMFTransform::MakeMatrixForReg( double *A, double *F, double *M, int *IM, int hn )

{

int i,j,k;

memset( A, 0, hn*(hn+1)/2*sizeof(double) );

memset( F, 0, hn*sizeof(double) );

for( k = 0; k < hn; k++ )

for( j = 0; j < IM[k+1]-IM[k]; j++ )

for( i = 0; i <= j; i++ )

A[(j+k)*(j+k+1)/2+i+k] += M[IM[k]+i]*M[IM[k]+j];

for( k = 0; k < hn; k++ )

for( j = 0; j < IM[k+1]-IM[k]; j++ )

F[j+k] += M[IM[k]+j]*EMF[k];

}

int CEMFTransform::RestHvsLin( double *HT, double *HEMF, int &hn )

{

int i, j;

// Memory map

//  int IM[n]              Profile for M

//  double M[n*(n+1)/2]      SLAE matrix

// Total memory n*(n+1)/2*sizeof(double)+n*sizeof(int)

double *M;

int *IM;

int res;

res = DelImpulsTime();

if( res != errOk )return res;

if( m < 2 || n < 2 )return errTooFewPoints;

IM = (int *)&IntMem[2*m*sizeof(double)];

M = (double *)&IntMem[2*m*sizeof(double)+(n+2)*sizeof(int)];

res = MakeMatrix( M, IM, HT, hn );

if( res != errOk )return res;

// Solve SLAE

double sum, maxelem = fabs(M[0]);

/* debug information */

/*

ofstream dofs( "a.dat" );

dofs.precision( 3 );

dofs.setf( ios::scientific );

for( i = 0; i < hn; i++ )

{

for( j = 0; j < IM[i+1]-IM[i]; j++ )dofs << M[IM[i]+j] << ' ';

dofs << endl;

}

dofs << endl;

dofs.precision( 6 );

dofs << "DE:\tmaxNDE:\tPos:\tNDE/DE:" << endl;

for( i = 0; i < hn-1; i++ )

{

maxelem = fabs(M[IM[i]+1]);

res = i+1;

for( j = 2; j < IM[i+1]-IM[i]; j++ )

if( maxelem < fabs(M[IM[i]+j]) )

{

maxelem = fabs(M[IM[i]+j]);

res = j+i;

}

dofs << M[IM[i]] << '\t' << maxelem << '\t' << res << '/' << (i+IM[i+1]-1-IM[i]) << '\t' << (maxelem/M[IM[i]]) << endl;

}

dofs.close();*/

// -------------------for( i = 1; i <= IM[hn-1]; i++ )

maxelem = max( fabs(M[i]), maxelem );

maxelem *= 1e-15;

for( i = hn-1; i >= 0; i-- )