Skip to content
Snippets Groups Projects
Commit dc507bcf authored by Wuttke, Joachim's avatar Wuttke, Joachim
Browse files

corr & debmsg

parent b23a87f1
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
pub/src/plot.cpp
+ 83
− 79
View file @ dc507bcf
......@@ -53,23 +53,23 @@ void NPlot::Plot( class CPlot *plot, bool add )
POld fd;
while((fd=fiter())) {
JSel.evaluate( 0, fd->nJ()-1 );
for( uint iv=0; iv<JSel.size(); ++iv ) {
for ( uint iv=0; iv<JSel.size(); ++iv ) {
vx = &( fd->VS( JSel.V[iv] )->x );
for( uint i=0; i<vx->size(); ++i ){
if( plot->X.logflag && (*vx)[i]<=0 )
for ( uint i=0; i<vx->size(); ++i ){
if ( plot->X.logflag && (*vx)[i]<=0 )
continue;
if( first ){
if ( first ){
inf = (*vx)[i];
sup = (*vx)[i];
first = false;
} else {
if( (*vx)[i]<inf ) inf = (*vx)[i];
if( (*vx)[i]>sup ) sup = (*vx)[i];
if ( (*vx)[i]<inf ) inf = (*vx)[i];
if ( (*vx)[i]>sup ) sup = (*vx)[i];
}
}
}
}
if( first )
if ( first )
throw string( "x range is empty" );
plot->X.setRoundedLimits( "x", inf, sup );
}
......@@ -83,13 +83,13 @@ void NPlot::Plot( class CPlot *plot, bool add )
while( (f=iter()) ){
uint k = iter.k();
POlc fc = P2C( f );
if( fc && fc->kconv!=-1 )
if ( fc && fc->kconv!=-1 )
continue;
POld fd = P2D( f );
JSel.evaluate( 0, f->nJ()-1 );
for( uint iv=0; iv<JSel.size(); ++iv ){
for ( uint iv=0; iv<JSel.size(); ++iv ){
uint j = JSel.V[iv];
if( fd ){
if ( fd ){
s = fd->VS(j);
} else {
S->x.resize( npts );
......@@ -99,9 +99,9 @@ void NPlot::Plot( class CPlot *plot, bool add )
s = S;
}
for ( uint i=0; i<s->size(); i++ ) {
if( !plot->X.contains( s->x[i] ) )
if ( !plot->X.contains( s->x[i] ) )
continue;
if( plot->Y.logflag && s->y[i]<=0 )
if ( plot->Y.logflag && s->y[i]<=0 )
continue;
if (first) {
inf=s->y[i], sup=s->y[i];
......@@ -113,7 +113,7 @@ void NPlot::Plot( class CPlot *plot, bool add )
}
}
}
if( first )
if ( first )
throw string( "y range is empty" );
plot->Y.setRoundedLimits( "y", inf, sup );
}
......@@ -134,7 +134,6 @@ void NPlot::Plot( class CPlot *plot, bool add )
uint k = fiter.k();
POld fd = P2D( f );
POlc fc = P2C( f );
int npts = 0;
JSel.evaluate( 0, f->nJ()-1 );
for ( uint iv=0; iv<JSel.size(); ++iv){
......@@ -151,7 +150,7 @@ void NPlot::Plot( class CPlot *plot, bool add )
vector<double> xp, yp, dyp;
for ( uint i=0; i<n; i++ ) {
double x = s->x[i];
if( !plot->X.contains(x) ) {
if ( !plot->X.contains(x) ) {
if ( x<=plot->X.inf ){
x = plot->X.inf;
nxl++;
......@@ -164,7 +163,7 @@ void NPlot::Plot( class CPlot *plot, bool add )
continue;
}
double y = s->y[i];
if( !plot->Y.contains(y) ) {
if ( !plot->Y.contains(y) ) {
if ( y<=plot->Y.inf ){
y = plot->Y.inf;
nyl++;
......@@ -176,35 +175,34 @@ void NPlot::Plot( class CPlot *plot, bool add )
if ( !plot->Y.force )
continue;
}
if( np==plot->maxpoints && np<n ) {
if ( np==plot->maxpoints && np<n ) {
cout << "reached maxpoints at " << s->x[i] << "\n";
}
xp.push_back( x );
yp.push_back( y );
if( s->dy.size() )
if ( s->dy.size() )
dyp.push_back( s->dy[i] );
np = xp.size();
if( np > plot->maxpoints ) // continue in large steps
if ( np > plot->maxpoints ) // continue in large steps
i += plot->maxpoints;
}
if ( np==0 ) {
cout << "file " << k << " spec " << j << " all " << n <<
" points outside plot range:\n";
if( nxl )
if ( nxl )
cout << " " << nxl << " points < xmin\n";
if( nxh )
if ( nxh )
cout << " " << nxh << " points < xmax\n";
if( nyl )
if ( nyl )
cout << " " << nyl << " points < ymin\n";
if( nyh )
if ( nyh )
cout << " " << nyh << " points > ymax\n";
continue;
} else {
plot->addSpec( false, pstyle++,
xp, yp, dyp, f->V[j]->z,
f->xco.str(), f->yco.str(),
"data file " + strg(k) +
" spectrum "+strg(j) );
plot->addSpec(
false, pstyle++, xp, yp, dyp, f->V[j]->z,
f->xco.str(), f->yco.str(),
"data file " + strg(k) + " spectrum "+strg(j) );
}
} else if ( fc && fc->kconv!=-1 ) { // plot curve
......@@ -216,21 +214,21 @@ void NPlot::Plot( class CPlot *plot, bool add )
if ( !fconv )
throw string( "convolutand is not a data file" );
vector<double> xp, yp;
for( uint i=0; i<fconv->VS(jconv)->size(); ++i ){
for ( uint i=0; i<fconv->VS(jconv)->size(); ++i ){
double x = fconv->VS(jconv)->x[i];
if( plot->X.contains( x ) )
if ( plot->X.contains( x ) )
xp.push_back( x );
}
fc->curve_val_vec( &yp, xp, k, j );
vector<double> xo, yo;
for( uint i=0; i<xp.size(); ++i ){
if( !( plot->X.contains( xp[i] ) &&
for ( uint i=0; i<xp.size(); ++i ){
if ( !( plot->X.contains( xp[i] ) &&
plot->Y.contains( yp[i] ) ) )
continue;
xo.push_back( xp[i] );
yo.push_back( yp[i] );
}
if( xo.size()==0 ){
if ( xo.size()==0 ){
cout << "curve k="<<strg(k)<<", j="<<strg(j)<<" is empty\n";
return;
}
......@@ -243,11 +241,11 @@ void NPlot::Plot( class CPlot *plot, bool add )
uint npts = plot->equipoints;
vector<double> xc(npts), yc, xp, yp;
for( uint i=0; i<npts; ++i )
for ( uint i=0; i<npts; ++i )
xc[i] = plot->X.plotcoord2value( i/(npts-1.0) );
fc->curve_val_vec( &yc, xc, k, j );
for( uint i=0; i<npts; ++i ) {
if( plot->X.contains(xc[i]) && plot->Y.contains(yc[i]) ) {
for ( uint i=0; i<npts; ++i ) {
if ( plot->X.contains(xc[i]) && plot->Y.contains(yc[i]) ) {
xp.push_back( xc[i] );
yp.push_back( yc[i] );
}
......@@ -263,47 +261,35 @@ void NPlot::Plot( class CPlot *plot, bool add )
// - interpolate at discontinuities
uint npts = plot->equipoints;
double xhig = plot->X.plotcoord2value( 1 );
double xlow = plot->X.plotcoord2value( 0 );
bool logx = plot->X.logflag;
vector< vector<double> > xpv, ypv;
vector<double> xc, yc, xn, yn, xp, yp;
double dx = logx ?
exp( log( xhig/xlow ) / (npts-1.0) ) :
(xhig-xlow) / (npts-1.0);
bool plotfirst = true;
bool plotlast = false;
bool plotend = false;
bool started = false;
bool ended = false;
double xt, yt; //current point (i)
double xtm,ytm;// // point (i - 1)
// double xtmo,ytmo;// back up of (i - 1)
vector<double> xc, yc, xn, yn;
// start with equidistant grid:
for( uint i=0; i<npts; ++i )
for ( uint i=0; i<npts; ++i )
xc.push_back( plot->X.plotcoord2value( i/(npts-1.0) ) );
// refinement loop:
for( int iref=0; ; ++iref ) {
for ( int iref=0; ; ++iref ) {
cout << "DEB iref=" << iref << "\n";
// evaluate curve for grid xc:
fc->curve_val_vec( &yc, xc, k, j );
cout << "DEB evaluation -> " << xc.size() << "\n";
for( uint i=0; i<xc.size(); ++i )
for ( uint i=0; i<xc.size(); ++i )
printf ("DEB xy %12.6g %12.6g\n", xc[i], yc[i] );
// merge xc,yc and xn,yn:
if( iref==0 ) {
if ( iref==0 ) {
xn = xc; yn = yc;
} else {
vector<double> xa, ya;
uint ic=0, in=0;
while( ic<xc.size() || in<xn.size() ) {
if( in==xn.size() ||
ic!=xc.size() && xc[ic]<xn[in] ){
for ( ; ic<xc.size() && in<xn.size(); ) {
if ( fabs( xc[ic]-xn[ic] ) < 1e-7 ) {
++in; // skip new point
} else if ( xc[ic]<xn[in] ) {
xa.push_back( xc[ic] );
ya.push_back( yc[ic] );
++ic;
......@@ -313,27 +299,33 @@ void NPlot::Plot( class CPlot *plot, bool add )
++in;
}
}
for ( ; ic<xc.size(); ++ic ) {
xa.push_back( xc[ic] );
ya.push_back( yc[ic] );
}
for ( ; in<xn.size(); ++in ) {
xa.push_back( xn[in] );
ya.push_back( yn[in] );
}
xn = xa; yn = ya;
}
cout << "DEB merger -> " << xn.size() << "\n";
// end of loop
if( iref>=10 || xn.size()>=plot->maxpoints )
if ( iref>=10 || xn.size()>=plot->maxpoints )
break;
// eliminate points outside y range (except border points):
vector<double> xa, ya;
uint i=0;
while( i<npts ) {
while ( i<npts ) {
cout << "DEB eli i " << i << "\n";
xa.push_back( xn[i] );
ya.push_back( yn[i] );
if( !plot->Y.contains(yc[i]) ) {
for( uint ii=i+1; ii<npts; ++ii ) {
cout << "DEB eli ii " << ii << "\n";
if( plot->Y.contains(yc[ii]) ) {
if ( !plot->Y.contains(yc[i]) ) {
for ( uint ii=i+1; ii<npts; ++ii ) {
if ( plot->Y.contains(yc[ii]) ) {
i = ii>i+1 ? ii-1 : ii;
cout << "DEB eli new i " << i << "\n";
break;
}
}
......@@ -345,39 +337,52 @@ void NPlot::Plot( class CPlot *plot, bool add )
cout << "DEB within frame -> " << xn.size() << "\n";
// set additional base points:
bool insert_next = false;
xc.clear();
for( uint i=0; i<npts-1; ++i ) {
for ( uint i=0; i<npts-1; ++i ) {
if ( i<npts-2 ) {
double yi = ( (xn[i+1]-xn[i])*yn[i-1] +
(xn[i]-xn[i-1])*yn[i+1] ) /
(xn[i+1]-xn[i-1]);
if( fabs( yn[i+1] - yi ) > 0.01*
(plot->Y.sup - plot->Y.inf) ) {
xc.push_back( (xn[i-1] + 2*xn[i])/3 );
xc.push_back( (2*xn[i] + xn[i+1])/3 );
if ( fabs( yn[i] - yi ) > 0.01*
(plot->Y.sup - plot->Y.inf) ) {
xc.push_back( (xn[i-1] + xn[i])/2 );
cout << "DEB disc " << xc.last() << "\n";
insert_next = true;
continue;
}
}
if ( insert_next ) {
xc.push_back( (xn[i-1] + xn[i])/2 );
cout << "DEB next " << xc.last() << "\n";
insert_next = false;
continue;
}
if ( ( plot->Y.contains(yc[i]) &&
!plot->Y.contains(yc[i-1]) ) ||
( !plot->Y.contains(yc[i]) &&
plot->Y.contains(yc[i-1]) ) ) {
xc.push_back( (2*xn[i-1] + xn[i])/3 );
xc.push_back( (xn[i-1] + 2*xn[i])/3 );
cout << "DEB edge " << xc.last() << "\n";
xc.push_back( (xn[i-1] + xn[i])/2 );
continue;
}
}
// temporary check
for ( uint i=0; i<xc.size()-1; ++i )
if ( xc[i+1]<=xc[i] )
throw "BUG: new base points not sorted";
}
// divide into subcurves that do not leave the frame:
vector<double> xa, ya;
for( uint i=0; i<npts; ++i ) {
if( plot->Y.contains(yc[i]) ) {
for ( uint i=0; i<npts; ++i ) {
if ( plot->Y.contains(yc[i]) ) {
xa.push_back( xn[i] );
ya.push_back( yn[i] );
}
if( !plot->Y.contains(yc[i]) || i==npts-1 ) {
if( xa.size() ){
if ( !plot->Y.contains(yc[i]) || i==npts-1 ) {
if ( xa.size() ){
cout << "DEB segment -> " << xa.size() << "\n";
xpv.push_back(xa);
ypv.push_back(ya);
......@@ -390,11 +395,10 @@ void NPlot::Plot( class CPlot *plot, bool add )
// plot subcurves: TODO merge with above, eliminate xpv,ypv
for ( uint ipv=0; ipv<xpv.size();++ipv) {
if(!xpv[ipv].empty()){
if (!xpv[ipv].empty()) {
plot->addSpec(
true, cstyle, xpv[ipv], ypv[ipv],
novec, f->V[j]->z,
f->xco.str(), f->yco.str(),
novec, f->V[j]->z, f->xco.str(), f->yco.str(),
"curve file "+strg(k)+ " spectrum "+strg(j) );
}
}
......@@ -403,7 +407,7 @@ void NPlot::Plot( class CPlot *plot, bool add )
throw string( "PROGRAM ERROR plot invalid *e" );
}
plot->plotDoclines( f->lDoc );
if( fc )
if ( fc )
plot->plotDoclines( fc->pInfo() );
}
plot->showSpecs();
......
pub/test/test_curve_ci
+ 0
− 0
View file @ dc507bcf
No preview for this file type
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment