From a4de122efcdb24ffeaa5155fe29b3b0179a643a7 Mon Sep 17 00:00:00 2001
From: "Joachim Wuttke (laptop)" <j.wuttke@fz-juelich.de>
Date: Fri, 26 Mar 2010 10:04:46 +0100
Subject: [PATCH] cleanup (edif.cpp partly -> plot.cpp; plotaux.cpp merged
->dualplot.cpp)
---
pub/src/Makefile.am | 6 +-
pub/src/Makefile.in | 5 +-
pub/src/coord.cpp | 2 +-
pub/src/dualplot.cpp | 238 ++++++++++++++--------
pub/src/dualplot.h | 13 +-
pub/src/edif.cpp | 442 +----------------------------------------
pub/src/edif.h | 34 ++--
pub/src/frida2.cpp | 35 ++--
pub/src/manip.cpp | 1 -
pub/src/mystd.cpp | 8 -
pub/src/mystd.h | 5 -
pub/src/plot.cpp | 463 +++++++++++++++++++++++++++++++++++++++++++
pub/src/plot.h | 3 +
pub/src/plotaux.cpp | 138 -------------
pub/src/plotaux.h | 6 -
pub/src/rssm.cpp | 1 +
16 files changed, 666 insertions(+), 734 deletions(-)
create mode 100644 pub/src/plot.cpp
create mode 100644 pub/src/plot.h
delete mode 100644 pub/src/plotaux.cpp
delete mode 100644 pub/src/plotaux.h
diff --git a/pub/src/Makefile.am b/pub/src/Makefile.am
index 5efbad83..8fce6fea 100644
--- a/pub/src/Makefile.am
+++ b/pub/src/Makefile.am
@@ -44,8 +44,8 @@ olm.cpp \
olm.h \
opr.cpp \
opr.h \
-plotaux.cpp \
-plotaux.h \
+plot.cpp \
+plot.h \
rng.cpp \
rng.h \
rssm.cpp \
@@ -57,7 +57,7 @@ special.h \
xax_lex.h \
xax_lex.lpp \
xax_yacc.ypp
-AM_LDFLAGS = -lboost_regex -lread-plus -llmmin -lkww -lyaml-cpp
+AM_LDFLAGS = -lboost_regex -lread-plus -llmmin -lkww -lyaml-cpp
# frida2_LDADD $(BOOST_REGEX_LIB)
# -lyaml-cpp
diff --git a/pub/src/Makefile.in b/pub/src/Makefile.in
index ffdfcf5c..08957b30 100644
--- a/pub/src/Makefile.in
+++ b/pub/src/Makefile.in
@@ -57,7 +57,7 @@ am_frida2_OBJECTS = coord.$(OBJEXT) curve.$(OBJEXT) dualplot.$(OBJEXT) \
edif.$(OBJEXT) expr.$(OBJEXT) file_in.$(OBJEXT) \
file_out.$(OBJEXT) frida2.$(OBJEXT) func.$(OBJEXT) \
index.$(OBJEXT) list.$(OBJEXT) manip.$(OBJEXT) mystd.$(OBJEXT) \
- olm.$(OBJEXT) opr.$(OBJEXT) plotaux.$(OBJEXT) rng.$(OBJEXT) \
+ olm.$(OBJEXT) opr.$(OBJEXT) rng.$(OBJEXT) \
rssm.$(OBJEXT) spec.$(OBJEXT) special2.$(OBJEXT) \
xax_lex.$(OBJEXT) xax_yacc.$(OBJEXT)
frida2_OBJECTS = $(am_frida2_OBJECTS)
@@ -266,8 +266,6 @@ olm.cpp \
olm.h \
opr.cpp \
opr.h \
-plotaux.cpp \
-plotaux.h \
rng.cpp \
rng.h \
rssm.cpp \
@@ -389,7 +387,6 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/mystd.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/olm.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/opr.Po@am__quote@
-@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/plotaux.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/rng.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/rssm.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/spec.Po@am__quote@
diff --git a/pub/src/coord.cpp b/pub/src/coord.cpp
index c3552c5b..270867b3 100644
--- a/pub/src/coord.cpp
+++ b/pub/src/coord.cpp
@@ -45,7 +45,7 @@ char CoordDecode( const string& resp, CCoord *rval )
}
if (np!=0) return 'h';
- for (i=j; i>0 && strchr(" \t",resp[i-1]); i--);
+ for (i=j; i>0 && strchr(" \t",resp[i-1]); i--) ;
if (i==0) return 'h';
*rval = CCoord( resp.substr(0,i), resp.substr(j+1, n-j-2) );
diff --git a/pub/src/dualplot.cpp b/pub/src/dualplot.cpp
index 3af4d865..8b17548c 100644
--- a/pub/src/dualplot.cpp
+++ b/pub/src/dualplot.cpp
@@ -9,17 +9,10 @@
#include <string.h>
#include <iostream>
#include <math.h>
-#include <sys/types.h>
-#include <sys/stat.h>
-#include <unistd.h>
#include <fcntl.h>
#include <boost/format.hpp>
#include "mystd.h"
-#include "plotaux.h"
-#include "coord.h"
-#include <ask_simple_value.h>
-#include <ask_and_callback.h>
#include "dualplot.h"
using boost::format;
@@ -73,6 +66,17 @@ void CAxis::set_from_string( string in )
throw string( "invalid lower bound, expecting real number" );
if ( !mystd::any2dbl(s2,&sup) )
throw string( "invalid upper bound, expecting real number" );
+ // check new range
+ if ( logflag && inf<=0 && inf!=-INFINITY ) {
+ setAuto();
+ throw string( "log scale requires range above 0" );
+ }
+ if (!logflag && ( ( inf!=-INFINITY && inf<-3e38 ) ||
+ ( sup!=+INFINITY && sup>+3e38 ) ) ) {
+ setAuto();
+ throw string( "lin scale exceeds typical PS implementation range "
+ "(abs<3e38)" );
+ }
}
//! Replace limits by rounded values.
@@ -182,31 +186,6 @@ string CAxis::str() const
return strg(inf) + " " + strg(sup);
}
-//! Query and check range.
-
-void CAxis::Ask( string quest )
-{
- string def, in;
- // query input
- def = str();
- in = sask( quest, def );
- if (in==def)
- return;
- if (in!="")
- set_from_string( in );
- // check new range
- if (logflag && inf<=0 && inf!=-INFINITY) {
- setAuto();
- throw string( "log scale requires range above 0" );
- }
- if (!logflag && ( ( inf!=-INFINITY && inf<-3e38 ) ||
- ( sup!=+INFINITY && sup>+3e38 ) ) ) {
- setAuto();
- throw string( "lin scale exceeds typical PS implementation range "
- "(abs<3e38)" );
- }
-}
-
//! Describe range, for use in table of plot frames.
string CAxis::info() const
@@ -225,6 +204,128 @@ double CAxis::pc( double v ) const
return CPLOT_PSMAX * value2plotcoord( v );
}
+void CAxis::calc_ticks( int *ntack, double **tack, int *ntpt, double *ticklim)
+ const
+{
+ int ir, nt, i;
+ double r, rd, rdr, d, vsub, vsup, tack0, tackn, dtack;
+
+ if ( inf >= sup )
+ throw "PROGRAM ERROR/ Ticks/ Bad Limits: " +
+ strg(inf) + ", " + strg(sup);
+
+ if (!logflag) {
+ r = log10(sup - inf);
+ ir = (int) ((r<0) ? r-1 : r);
+ rd = r - ir; // fractional part of r
+ rdr = pow(10., rd);
+ if (rdr > 9.99999) {
+ dtack = 2.5 * pow(10., ir); // spacing between tacks
+ *ntpt = 5; // ticks per tack
+ } else if (rdr > 5.00001) {
+ dtack = 2 * pow(10., ir);
+ *ntpt = 4;
+ } else if (rdr > 2.500005) {
+ dtack = 1 * pow(10., ir);
+ *ntpt = 5;
+ } else if (rdr > 1.600003) {
+ dtack = .5 * pow(10., ir);
+ *ntpt = 5;
+ } else if (rdr > 1.250002) {
+ dtack = .4 * pow(10., ir);
+ *ntpt = 4;
+ } else {
+ dtack = .25 * pow(10., ir);
+ *ntpt = 5;
+ }
+ d = inf / dtack;
+ tack0 = dtack * (int) (d<=0?d-1:d);
+ d = sup / dtack;
+ tackn = dtack * ((int) (d<0?d-1:d) + 1);
+ vsub = inf - (sup-inf)*1.e-4;
+ vsup = sup + (sup-inf)*1.e-4;
+ nt = (int) ((tackn-tack0) / dtack) + 3; // allocate enough
+ *tack = (double*) malloc(sizeof(double)*nt);
+ *ntack = 0;
+ for (i=0; i<nt; i++) {
+ d = tack0 + i*dtack;
+ if (vsub<=d && d<=vsup)
+ (*tack)[(*ntack)++] = d;
+ }
+
+ ticklim[0] = tack0;
+ ticklim[1] = (*tack)[(*ntack)-1] + dtack;
+
+ } else { // log scale
+
+ static double eins = 1.0000009999999999;
+ double rlgmin, rlgmax, rlgrel, d1, d2, r0;
+ int minexp, maxexp, incr;
+
+ if (inf <= 0)
+ throw string( "PROGRAM ERROR/ Ticks/ negative log argument" );
+
+ rlgmin = log10(inf);
+ rlgmax = log10(sup);
+ rlgrel = rlgmax - rlgmin;
+ if (rlgrel > 40)
+ throw string( "Excessive log range" );
+
+ minexp = (rlgmin>0 ?
+ (int) (rlgmin+1e-6) :
+ ((int) (rlgmin-1e-6) -1));
+ maxexp = (rlgmax>0 ?
+ (int) (rlgmax+1e-6) :
+ ((int) (rlgmax-1e-6) -1));
+
+ incr = (int) (rlgrel / 7 + 1); // #decades per tack
+ if (incr==2) incr = 1; // preferences: 10^1, 10^3, 10^6 etc.
+ if (incr>3) incr = 3 * (incr % 3);
+
+ if (incr > 1) {
+ minexp -= minexp % incr;
+ maxexp -= maxexp % incr;
+ }
+
+ /* - Set large ticks array : */
+ static int mt=30;
+ *tack = (double*) malloc(sizeof(double)*mt);
+
+ *ntack = 0;
+ d1 = inf / eins;
+ d2 = sup * eins;
+ for ( i = minexp; incr < 0 ? i >= maxexp : i <= maxexp; i += incr ) {
+ r0 = pow(10., i);
+ if (d1<=r0 && r0<=d2) {
+ if(*ntack>=mt)
+ throw string( "PROG ERR too many tacks" );
+ (*tack)[(*ntack)++] = r0;
+ }
+ }
+
+ if (*ntack >= 1) {
+ ticklim[0] = (*tack)[0] / pow(10., incr);
+ ticklim[1] = (*tack)[(*ntack)-1] * pow(10., incr);
+ } else {
+ // untested. just to allow plots when no tacks in window
+ ticklim[0] = pow( 10., minexp );
+ ticklim[1] = pow( 10., maxexp+1 );
+ }
+
+ /** determine # ticks / tack **/
+ if (incr==1 && rlgrel<=7) {
+ *ntpt = 9;
+ } else if (incr==1 && rlgrel<=7) {
+ *ntpt = 3; // 1-2-5-10 - Schritte
+ } else if (incr==1) {
+ *ntpt = 1;
+ } else {
+ *ntpt = -incr;
+ // means incr ticks per tack
+ // means 1 tick per decade
+ }
+ }
+}
//**************************************************************************//
//* CPlot *//
@@ -241,22 +342,18 @@ CPlot::CPlot( uint _iPlot, bool _logx, bool _logy ) :
string fn_fifo = string("/tmp/gnuplot-") + getenv( "LOGNAME" );
mystd::system( "rm -f "+fn_fifo );
- if (mkfifo( fn_fifo.c_str(), 0666 )) {
- printf("SEVERE SYSTEM ERROR cannot make fifo %s"
- " - will not be able to print\n", fn_fifo.c_str() );
- return;
- }
+ if (mkfifo( fn_fifo.c_str(), 0666 ))
+ throw "SEVERE SYSTEM ERROR cannot make fifo " + fn_fifo +
+ ": will not be able to print";
mystd::system( "gnuplot -title " + strg(iPlot) + " -noraise" +
" < " + fn_fifo + " &" );
// we use open instead of fopen or ofstream,
// because we need non-blocking mode.
- if (!(gp_fifo = open(fn_fifo.c_str(), O_WRONLY))) {
- printf("SEVERE SYSTEM ERROR cannot open gp_fifo"
- " - will not be able to print\n");
- return;
- }
+ if (!(gp_fifo = open(fn_fifo.c_str(), O_WRONLY)))
+ throw "SEVERE SYSTEM ERROR cannot open fifo " + fn_fifo +
+ ": will not be able to print";
fcntl(gp_fifo,F_SETFL,O_NONBLOCK);
gp_write( string("set terminal x11") );
@@ -286,9 +383,7 @@ void CPlot::clearFrame()
void CPlot::setAux( string cmd )
{
- if ( cmd=="gnp" )
- maxpoints = iask("Max # points in plot", maxpoints);
- else if ( cmd=="gxl" ) {
+ if ( cmd=="gxl" ) {
X.setLog( !X.logflag );
printf( "set x %s\n", X.logflag ? "log" : "lin" ); }
else if ( cmd=="gxl+" )
@@ -322,7 +417,7 @@ void CPlot::setAux( string cmd )
//! Plot coordinate frame (axes, ticks, labels).
-void CPlot::plotFrame()
+void CPlot::plotFrame( string xlabel, string ylabel )
{
gp_write( "set nologscale" );
string whichlog="";
@@ -348,11 +443,10 @@ void CPlot::plotFrame()
ps_accu.push_back( "\n/xPlotFrame {" );
if ( X.logflag && X.inf<= 0 )
throw "BUG: x log incompatible with limits " + X.str();
- plotaux::calc_ticks( X.logflag, X.inf, X.sup,
- &ntack, &tack, &ntpt, ticklim );
+ X.calc_ticks( &ntack, &tack, &ntpt, ticklim );
ps_ticktack(ntack, tack, ntpt, ticklim, &X);
free(tack);
- snprintf( outlin, CPLOT_LINSIZ-4, " {(%s", X.C.ps_str().c_str() );
+ snprintf( outlin, CPLOT_LINSIZ-4, " {(%s", xlabel.c_str() );
strncat( outlin, ")}\n", CPLOT_LINSIZ );
ps_accu.push_back( outlin );
ps_accu.push_back( " 0 10 0 0 0 90 "
@@ -365,11 +459,10 @@ void CPlot::plotFrame()
ps_accu.push_back( "\n/yPlotFrame {" );
if ( Y.logflag && Y.inf<= 0 )
throw "BUG: y log incompatible with limits " + Y.str();
- plotaux::calc_ticks( Y.logflag, Y.inf, Y.sup,
- &ntack, &tack, &ntpt, ticklim );
+ Y.calc_ticks( &ntack, &tack, &ntpt, ticklim );
ps_ticktack(ntack, tack, ntpt, ticklim, &Y);
free(tack);
- snprintf( outlin, CPLOT_LINSIZ, " {(%s", Y.C.ps_str().c_str() );
+ snprintf( outlin, CPLOT_LINSIZ, " {(%s", ylabel.c_str() );
strncat( outlin, ")}\n", CPLOT_LINSIZ );
ps_accu.push_back( outlin );
ps_accu.push_back( " 0 10 0 0 90 0 "
@@ -481,16 +574,8 @@ void CPlot::plotDoclines( const vector<string>& lDoc )
//! Write buffered plot to postscript file.
-void CPlot::writePostscript( bool full_outfile )
+void CPlot::writePostscript( string ps_outdir, string ps_head, string ps_dict )
{
- string ps_outdir, ps_head, ps_dict;
-
- // read configuration parameters:
- ps_outdir = NRead::wofmac( "\\psdir" );
- ps_head = NRead::wofmac( "\\pshead" );
- if ( full_outfile )
- ps_dict = NRead::wofmac( "\\psdict" );
-
// construct output file name:
FILE *pssav;
string flong, cmd, outf;
@@ -498,7 +583,7 @@ void CPlot::writePostscript( bool full_outfile )
if (ps_fnum>=999)
throw string( "graph file number overflow" );
outf = ps_outdir + str( format( "l%d." ) % ++ps_fnum ) +
- ( full_outfile ? "ps" : "psa" );
+ ( ps_dict=="" ? "psa" : "ps" );
if (!(pssav = mystd::glob_fopen(outf.c_str(), "", "", "r")))
break; // legal exit
fclose(pssav);
@@ -506,22 +591,20 @@ void CPlot::writePostscript( bool full_outfile )
cout << "save plot in " << outf << "\n";
// copy headers to output file:
- cmd = string("cat ") + ( full_outfile ? ps_dict : "" ) + " " +
- ps_head + " > " + outf + "\n";
+ cmd = string("cat ") +
+ ps_dict + " " + // ps_dict may be ""
+ ps_head + " > " +
+ outf + "\n";
mystd::system( cmd );
// a redundant check of existence of the outfile:
- if (!(pssav = mystd::glob_fopen(outf.c_str(), "", "", "r", &flong))) {
+ if (!(pssav = mystd::glob_fopen(outf.c_str(), "", "", "r", &flong)))
throw "have not created file " + outf;
- return;
- }
fclose(pssav);
// append specific output to output file:
- if (!(pssav = fopen(flong.c_str(), "a+"))) {
- printf ("cannot write (append) to file %s\n", flong.c_str());
- return;
- }
+ if (!(pssav = fopen(flong.c_str(), "a+")))
+ throw "cannot write (append) to file " + flong;
for( uint i=0; i<ps_accu.size(); ++i ){
// fprintf does not work here because output line may contain "%"
fwrite( ps_accu[i].c_str(), 1, ps_accu[i].size(), pssav );
@@ -541,19 +624,6 @@ void CPlot::writePostscript( bool full_outfile )
fclose(pssav);
}
-//! Endless loop: query for gnuplot commands.
-
-void CPlot::feedGnuplot()
-{
- string cmd;
- cout << "entering mygnuplot - to leave, type q\n";
- while(1) {
- cmd = sask("mygnuplot> ");
- if (cmd.substr(0,1)=="q") return;
- gp_write(cmd);
- }
-}
-
//! Info line to characterize this plot window.
string CPlot::info() const
diff --git a/pub/src/dualplot.h b/pub/src/dualplot.h
index 6654c636..339a18c4 100644
--- a/pub/src/dualplot.h
+++ b/pub/src/dualplot.h
@@ -1,8 +1,7 @@
class CAxis {
public:
double inf, sup;
- CCoord C;
- bool logflag; // 'log' is reserved for math.h
+ bool logflag;
bool force;
CAxis( bool _log ) :
@@ -24,9 +23,10 @@ class CAxis {
double value2plotcoord( double v ) const;
double plotcoord2value( double c ) const;
double inf_pos() const;
+ void calc_ticks( int *ntack, double **tack, int *ntpt, double *ticklim )
+ const;
};
-
#define CPLOT_PSMAX 10
#define CPLOT_LINSIZ 80
@@ -39,7 +39,7 @@ class CPlot {
CPlot( uint _iPlot, bool _logx, bool _logy );
void clearFrame();
- void plotFrame();
+ void plotFrame( string xlabel, string ylabel );
void addSpec( bool as_line, int style_no,
const vector<double>& xp,
const vector<double>& yp, const vector<double>& dyp,
@@ -47,12 +47,11 @@ class CPlot {
string xco, string yco, string info );
void showSpecs();
void plotDoclines( const vector<string>& lDoc );
- void writePostscript( bool full_outfile );
+ void writePostscript( string ps_outdir, string ps_head, string ps_dict );
void setAux( string cmd );
- void feedGnuplot();
string info() const;
-private: // can we move this to the .cpp file ?
+private:
int gp_fifo;
void gp_write( string );
int gp_fno;
diff --git a/pub/src/edif.cpp b/pub/src/edif.cpp
index 080d0333..8c418c42 100644
--- a/pub/src/edif.cpp
+++ b/pub/src/edif.cpp
@@ -1,6 +1,6 @@
//**************************************************************************//
//* FRIDA: fast reliable interactive data analysis *//
-//* edif.cpp: edit on-line files (text_edit, dir, make, plot, ...) *//
+//* edif.cpp: edit on-line files (text_edit, dir, make, ...) *//
//* (C) Joachim Wuttke 1990-, v2(C++) 2001- *//
//* http://frida.sourceforge.net *//
//**************************************************************************//
@@ -10,13 +10,10 @@
#include <iostream>
#include <algorithm>
-using namespace std;
-
#include "mystd.h"
#include "olm.h"
#include "expr.h"
#include <ask_simple_value.h>
-#include "dualplot.h"
#include "edif.h"
#include "index.h"
@@ -263,7 +260,6 @@ user_choice:
} else if (choice=="" || choice=="q") {
return;
} else if (Choice.parse(choice, 0, AllR.size())) {
- cout << mystd::BEL;
goto user_choice;
}
@@ -312,7 +308,6 @@ user_choice:
CParam(AllR[ii], val));
}
} else {
- cout << mystd::BEL;
goto action_on_par;
} // action on rpar
} // list of rpars
@@ -950,438 +945,3 @@ void NEdif::MakeGrid(void)
NOlm::SelNew();
NOlm::OloAdd(&olf);
}
-
-//**************************************************************************//
-//* Plotting *//
-//**************************************************************************//
-
-
-
-double NEdif::FindAngleCos(double x1,double y1, //point 1
- double x2,double y2, //point 2
- double x3,double y3 )//point 3
-// Li Huang mar10
-{
- /*
- (1) ________(2)
- \)cos
- \
- \(3)
-
- (x1,y1), (x2,y2)
- (x2,y2) and (x3,y3)
- a = x2 - x1, b = y2 - y1
- c = x3 - x2, d = y3 - y2.
- The so called dot product yields the equation
- Cosine( Theta ) = ( (a x c) + (b x d) )/(Sqrt(a^2 + b^2) x Sqrt(c^2 + d^2) ) = Z.
- */
- double a,b,c,d,cos;
- a = x2 - x1;
- b = y2 - y1;
- c = x3 - x2;
- d = y3 - y2;
- cos = ( (a * c) + (b * d) )/(sqrt(a*a + b*b) * sqrt(c*c + d*d) );
- return cos;
-}
-
-//! The next point to be plotted is returned as *xt, *yt.
-
-void NEdif::FindNextPoint( const COlc *fc,
- double xtm, double ytm, double dx,
- double* xt, double* yt, uint k, uint j )
-// Li Huang mar10
-{
- uint i = 0;
- double cos123 = 0;
- double cos234 = 0;
- double x1 = xtm;
- double y1 = ytm;
- double x2 = *xt;
- double y2 = *yt;
- double x3,y3;
- double x4,y4;
- /*
- /)o1__ (3)
- (2)_____________
- / \)o2
- / \
- (1)/ \(4)
-
- (xt,yt)__________(xt+dx,yt+dx)
- / \
- / \
- (xtm,ytm)/ \(xt+2*dx,yt+2*dx)
-
- o1 < 1/4 degree && o2 < 1/4 degree
- cos(1/4 degree) = 0.99999
- */
- do{
- x3 = x2 + dx;
- y3 = fc->T->tree_curve_val_sca( x3, k, j );
- x4 = x3 + dx;
- y4 = fc->T->tree_curve_val_sca( x4, k, j );
- cos123 = FindAngleCos(x1,y1,x2,y2,x3,y3);
- cos234 = FindAngleCos(x2,y2,x3,y3,x4,y4);
- dx = dx/2;
- i++;
- }while ( ( cos123 < 0.99999 || cos234 < 0.99999 ) && i <= 7);
- *xt = x3;
- *yt = y3;
-}
-
-void NEdif::InterY(double xtm,double ytm,double xt,double yt,
- double* xtt,double *ytt,double border)
-{
- *ytt = border;
- *xtt = (*ytt-ytm)*(xt-xtm)/(yt-ytm)+xtm;
-}
-
-
-//! Plot spectrum.
-
-void NEdif::Plot( CPlot *plot, bool add )
-{
- static CList JSel;
- CEle *e;
- uint k, j;
- const COld *fd;
- const COlc *fc;
- CSpec S;
- vector<double> novec;
- static int pstyle = 0, cstyle = 0;
- static string jSel = "";
- NOlm::JSelAsk( "Plot which spectra", &jSel, &JSel );
-
- if (!add) {
-
- pstyle = 0;
- cstyle = 0;
-
- // determine x-range ?
-
- double inf, sup;
- if (!plot->X.finite()) {
- NOlm::IterateD fiter;
- const vector<double> *vx;
- bool first = true;
- while((fd=fiter())) {
- JSel.evaluate( 0, fd->nSpec()-1 );
- 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 )
- continue;
- 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( first )
- throw string( "x range is empty" );
- plot->X.setRoundedLimits( "x", inf, sup );
- }
-
- // determine y-range ?
- if (!plot->Y.finite()) {
- bool first=true;
- NOlm::IterateEle eiter;
- int npts = 100;
- const CSpec *s;
- while( (e=eiter()) ){
- k = eiter.SelNo();
- fc = e->C();
- if( fc && fc->kconv!=-1 )
- continue;
- fd = e->D();
- JSel.evaluate( 0, e->nSpec()-1 );
- for( uint iv=0; iv<JSel.size(); ++iv ){
- uint j = JSel.V[iv];
- if( fd ){
- s = &(fd->VS[j]);
- } else {
- S.x.resize( npts );
- for ( uint i=0; i<npts; i++ )
- S.x[i] = plot->X.plotcoord2value( i/(npts-1.0) );
- fc->T->tree_curve_val_vec( &(S.y), S.x, k, j );
- s = &S;
- }
- for (uint i=0; i<s->size(); i++) {
- if( !plot->X.contains( s->x[i] ) )
- continue;
- if( plot->Y.logflag && s->y[i]<=0 )
- continue;
- if (first) {
- inf=s->y[i], sup=s->y[i];
- first = false;
- } else {
- if (s->y[i]<inf) inf=s->y[i];
- if (s->y[i]>sup) sup=s->y[i];
- }
- }
- }
- }
- if( first )
- throw string( "y range is empty" );
- plot->Y.setRoundedLimits( "y", inf, sup );
- }
-
- // build label: (zur Zeit nur vom ersten File; definiere +=)
- plot->X.C = NOlm::MOM[*(NOlm::FSel.V.begin())].P()->xco;
- plot->Y.C = NOlm::MOM[*(NOlm::FSel.V.begin())].P()->yco;
-
- // draw new frame:
- plot->clearFrame();
- plot->plotFrame();
- }
-
- // plot data:
- NOlm::IterateEle eiter;
- while( (e=eiter()) ){
- k = eiter.SelNo();
- fc = e->C();
- fd = e->D();
- JSel.evaluate( 0, e->nSpec()-1 );
- for ( uint iv=0; iv<JSel.size(); ++iv){
- j = JSel.V[iv];
- // cout << "DEBUG edif plot k="<<k<<", j="<<j<<"\n";
- if ( fd ) {
- // plot data
- const CSpec *s;
- s = &( fd->VS[j] );
- uint n=s->x.size();
- if ( n!=s->y.size() )
- throw string( "BUG: plot: x.size<>y.size" );
- // Filter: x,y -> xp,yp if inside frame
- uint np=0, nxl=0, nxh=0, nyl=0, nyh=0;
- vector<double> xp, yp, dyp;
- for ( uint i=0; i<n; i++ ) {
- double x = s->x[i];
- if( !plot->X.contains(x) ) {
- if ( x<=plot->X.inf ){
- x = plot->X.inf;
- nxl++;
- }
- if ( x>=plot->X.sup ){
- x = plot->X.sup;
- nxh++;
- }
- if ( !plot->X.force )
- continue;
- }
- double y = s->y[i];
- if( !plot->Y.contains(y) ) {
- if ( y<=plot->Y.inf ){
- y = plot->Y.inf;
- nyl++;
- }
- if ( y>=plot->Y.sup ){
- y = plot->Y.sup;
- nyh++;
- }
- if ( !plot->Y.force )
- continue;
- }
- if( np==plot->maxpoints && np<n ) {
- printf("reached maxpoints at %g\n", s->x[i]);
- }
- xp.push_back( x );
- yp.push_back( y );
- if( s->dy.size() )
- dyp.push_back( s->dy[i] );
- np = xp.size();
- if( np > plot->maxpoints ) // continue in large steps
- i += plot->maxpoints;
- }
- if ( np==0 ) {
- printf( "file %u spec %u:"
- " all %u points outside plot range:\n",
- k, j, n );
- if( nxl )
- printf( " %u points < xmin\n", nxl );
- if( nxh )
- printf( " %u points > xmax\n", nxh );
- if( nyl )
- printf( " %u points < ymin\n", nyl );
- if( nyh )
- printf( " %u points > ymax\n", nyh );
- continue;
- } else {
- plot->addSpec( false, pstyle++,
- xp, yp, dyp, *(e->Z(j)),
- e->P()->xco.str(), e->P()->yco.str(),
- "data file " + strg(k) + " spectrum "+strg(j) );
- }
- } else if ( fc && fc->kconv!=-1 ) {
- // plot curve, but only at x points of convolutand
- uint kconv = fc->kconv;
- uint jconv = CTree::setConvJ( k, j, kconv );
- COld *fconv = NOlm::MOM[kconv].D();
- if ( !fconv )
- throw string( "convolutand is not a data file" );
- vector<double> xp, yp;
- for( uint i=0; i<fconv->VS[jconv].size(); ++i ){
- double x = fconv->VS[jconv].x[i];
- if( plot->X.contains( x ) )
- xp.push_back( x );
- }
- fc->T->tree_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] ) &&
- plot->Y.contains( yp[i] ) ) )
- continue;
- xo.push_back( xp[i] );
- yo.push_back( yp[i] );
- }
- if( xo.size()==0 ){
- cout << "curve k="<<strg(k)<<", j="<<strg(j)<<" is empty\n";
- return;
- }
- plot->addSpec( true, cstyle++, xo, yo, novec, *(e->Z(j)),
- e->P()->xco.str(), e->P()->yco.str(),
- "curve file "+strg(k)+" spectrum "+strg(j) );
- } else if ( fc && false ) {
- // PRESENTLY UNUSED
- // plot ordinary curve
- // straight vectorial version, for rapid plotting
- int npts = 297;
- vector<double> xc(npts), yc, xp, yp;
- for( uint i=0; i<npts; ++i )
- xc[i] = plot->X.plotcoord2value( i/(npts-1.0) );
- fc->T->tree_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]) ) {
- xp.push_back( xc[i] );
- yp.push_back( yc[i] );
- }
- }
- plot->addSpec( true, cstyle++, xp, yp, novec, *(e->Z(j)),
- e->P()->xco.str(), e->P()->yco.str(),
- "curve file "+strg(k)+" spectrum "+strg(j) );
- } else if ( fc ) {
- // plot ordinary curve, taking care:
- // - stop/start line when leaving/reentering the frame (done)
- // - interpolate when crossing the horizontal border (done)
- // - interpolate at discontinuities (done)
- //iamhere
- uint npts = 97;
- double dx = plot->X.plotcoord2value( 1/(npts-1.0) ) -
- plot->X.plotcoord2value( 0/(npts-1.0) );
- vector< vector<double> > xpv, ypv;
- vector<double> xp, yp;
- 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)
- while ( (!plotend) ) {
- if( xpv.size()>180 )
- throw string( "curve oscillates too fast for "
- "interpolating plot mode" );
- if ( !plotlast ) {
- if ( plotfirst ){
- xt = plot->X.plotcoord2value( 0.0 );
- yt = fc->T->tree_curve_val_sca( xt, k, j );
- } else {
- FindNextPoint( fc, xtmo, ytmo, dx,
- &xt, &yt, k, j );
- }
- }
- if ( plot->X.contains(xt) ){// if xt <= X.sup
- if ( plot->Y.contains(yt) ) { //if in
- if ( plotfirst ) { //a new curve starts
- //with the first point as the starting point
- xp.push_back( xt );
- yp.push_back( yt );
- plotfirst = false;
- started = true;
- }
- else if ( !started ) { //a new curve starts
- //i.e. the case: out -> in
- double xtt,ytt;
- double border = ytm > plot->Y.sup ?
- plot->Y.sup : plot->Y.inf;
- //using the previous point and current point to
- //interpolate the intersection point
- InterY(xtm,ytm,xt,yt,&xtt,&ytt,border);
- xp.push_back( xtt );
- yp.push_back( ytt );
- xp.push_back( xt );
- yp.push_back( yt );
- started = true;
- }
- else if ( started ) { //continues a started curve
- //i.e. the case in -> in
- xp.push_back( xt );
- yp.push_back( yt );
- }
- }
- else if ( started ) {
- // a curve ends
- //i.e. the case: in -> in -> out
- double xtt,ytt;
- double border = yt > plot->Y.sup ?
- plot->Y.sup : plot->Y.inf;
- //using the previous point and current point to
- //interpolate the intersection point
- InterY(xtm,ytm,xt,yt,&xtt,&ytt,border);
- xp.push_back( xtt );
- yp.push_back( ytt );
- started = false;
- ended = true;
- }
- else if ( plotfirst ) { //a new curve starts
- //with the first point is not inside the plotting
- plotfirst = false;
- }
- if ( xt == plot->X.sup ) {
- ended = true;
- plotlast = true;
- plotend = true;
- }
- } else if ( ( xt > plot->X.sup ) && !plotlast ){
- xt = plot->X.sup;
- plotlast = true;
- } else
- throw string( "PROGRAM ERROR plot in fc: no ending" );
-
- xtmo = xtm;
- ytmo = ytm;
- xtm = xt;
- ytm = yt;
- if ( ended ){
- xpv.push_back(xp);
- ypv.push_back(yp);
- xp.clear();
- yp.clear();
- ended = false;
- }
- } // end loop over points
-
- for ( uint ipv=0; ipv<xpv.size();++ipv) {
- if(!xpv[ipv].empty()){
- plot->addSpec( true, cstyle, xpv[ipv], ypv[ipv],
- novec, *(e->Z(j)),
- e->P()->xco.str(), e->P()->yco.str(),
- "curve file "+strg(k)+" spectrum "+strg(j) );
- }
- }
- cstyle++;
- } else
- throw string( "PROGRAM ERROR plot invalid *e" );
- }
- plot->plotDoclines( e->P()->lDoc );
- if( fc )
- plot->plotDoclines( fc->pInfo() );
- }
- plot->showSpecs();
-}
diff --git a/pub/src/edif.h b/pub/src/edif.h
index a6c65fde..b045f22c 100644
--- a/pub/src/edif.h
+++ b/pub/src/edif.h
@@ -1,24 +1,14 @@
namespace NEdif {
- void Dir();
- void DirZ();
- void ListZ();
- void DirNumpar();
- void EditNumpar();
- void DirCoord();
- void EditFNam();
- void EditCoord( string which );
- void EditDoc();
- void ReadIn();
- void ReadTab( string qualif );
- void MakeGrid();
- void Plot( CPlot *plot, bool add );
-/* auxiliary for plotting: */
- double FindAngleCos(double x0,double y0,
- double x1,double y1,
- double x2,double y2 );
- void FindNextPoint( const COlc *fc,
- double xtm, double ytm, double dx,
- double* xt, double* yt, uint k, uint j );
- void InterY(double xtm,double ytm,double xt,double yt,
- double* xtt,double *ytt,double border);
+ void Dir();
+ void DirZ();
+ void ListZ();
+ void DirNumpar();
+ void EditNumpar();
+ void DirCoord();
+ void EditFNam();
+ void EditCoord( string which );
+ void EditDoc();
+ void ReadIn();
+ void ReadTab( string qualif );
+ void MakeGrid();
};
diff --git a/pub/src/frida2.cpp b/pub/src/frida2.cpp
index 3475cf19..442374e6 100644
--- a/pub/src/frida2.cpp
+++ b/pub/src/frida2.cpp
@@ -19,6 +19,7 @@
#include "dualplot.h"
#include "opr.h"
#include "edif.h"
+#include "plot.h"
#include "file_in.h"
#include "file_out.h"
#include "rssm.h"
@@ -283,16 +284,17 @@ int main()
" gf plot to short .psX file\n"
" gw list of plot windows\n"
" g<n> switch to plot window <n>\n"
- " gd dialog (within gnuplot)\n"
" p plot\n"
" pp plot with new autoranges\n"
" a add to plot\n"
;
} else if (cmd == "gx") {
- Plots[iPlot]->X.Ask("Set x range");
+ string range = sask( "Set x range" );
+ Plots[iPlot]->X.set_from_string( range );
} else if (cmd == "gy") {
- Plots[iPlot]->Y.Ask("Set y range");
+ string range = sask( "Set y range" );
+ Plots[iPlot]->Y.set_from_string( range );
} else if (cmd == "ga") {
Plots[iPlot]->X.setAuto();
Plots[iPlot]->Y.setAuto();
@@ -301,15 +303,20 @@ int main()
} else if (cmd == "gya") {
Plots[iPlot]->Y.setAuto();
} else if (cmd.substr(0,2) == "gx" ||
- cmd.substr(0,2) == "gy" ||
- cmd == "gnp" ) {
+ cmd.substr(0,2) == "gy" ) {
Plots[iPlot]->setAux( cmd );
- } else if (cmd == "gp") {
- Plots[iPlot]->writePostscript( true );
- } else if (cmd == "gf") {
- Plots[iPlot]->writePostscript( false );
- } else if (cmd == "gd") { // graph dialog
- Plots[iPlot]->feedGnuplot();
+ } else if ( cmd=="gnp" ) {
+ Plots[iPlot]->maxpoints =
+ iask( "Max # points in plot", Plots[iPlot]->maxpoints );
+ } else if (cmd == "gp" || cmd == "gf" ) {
+ string ps_outdir = NRead::wofmac( "\\psdir" );
+ string ps_head = NRead::wofmac( "\\pshead" );
+ string ps_dict;
+ if ( cmd=="gp" )
+ ps_dict = NRead::wofmac( "\\psdict" );
+ else
+ ps_dict = "";
+ Plots[iPlot]->writePostscript( ps_outdir, ps_head, ps_dict );
} else if (cmd == "gw") {
for( uint i=0; i<nPlot; ++i ){
printf( " %c %1i %s\n", i==iPlot ? '*' : ' ', i,
@@ -324,13 +331,13 @@ int main()
} else
throw string( "invalid graphic window number" );
} else if (cmd == "p") {
- NEdif::Plot( Plots[iPlot], 0);
+ NPlot::Plot( Plots[iPlot], 0);
} else if (cmd == "pp") {
Plots[iPlot]->X.setAuto();
Plots[iPlot]->Y.setAuto();
- NEdif::Plot( Plots[iPlot], 0 );
+ NPlot::Plot( Plots[iPlot], 0 );
} else if (cmd == "a") {
- NEdif::Plot( Plots[iPlot], 1 );
+ NPlot::Plot( Plots[iPlot], 1 );
} else if (cmd == "m") {
diff --git a/pub/src/manip.cpp b/pub/src/manip.cpp
index 72caa1ae..c4bdd286 100644
--- a/pub/src/manip.cpp
+++ b/pub/src/manip.cpp
@@ -18,7 +18,6 @@
#include "expr.h"
#include <ask_simple_value.h>
#include "manip.h"
-#include "dualplot.h" // wg. edif.h
#include "edif.h" // wg. ListZ()
#include "xax_lex.h"
#include "index.h"
diff --git a/pub/src/mystd.cpp b/pub/src/mystd.cpp
index c988eedf..8eee689e 100644
--- a/pub/src/mystd.cpp
+++ b/pub/src/mystd.cpp
@@ -24,14 +24,6 @@ namespace mystd_aux {
};
-//**************************************************************************//
-//* constants *//
-//**************************************************************************//
-
-const char mystd::BEL = 0x7;
-const double mystd::PI = 3.1415926535897931;
-
-
//**************************************************************************//
//* convert to string *//
//**************************************************************************//
diff --git a/pub/src/mystd.h b/pub/src/mystd.h
index 1439e600..da352ac0 100644
--- a/pub/src/mystd.h
+++ b/pub/src/mystd.h
@@ -1,5 +1,4 @@
#include <string>
-#include <fstream>
#include <vector>
#define SQR(a) ((a)*(a))
@@ -22,10 +21,6 @@ string yaml( const string );
namespace mystd {
- //! constants
- extern const double PI;
- extern const char BEL;
-
//! convert from string
bool str2vec( string inp, vector<double> *V, uint nmax=0, bool force=true );
bool any2dbl( const string& inp, double *val );
diff --git a/pub/src/plot.cpp b/pub/src/plot.cpp
new file mode 100644
index 00000000..39540bfa
--- /dev/null
+++ b/pub/src/plot.cpp
@@ -0,0 +1,463 @@
+//**************************************************************************//
+//* FRIDA: fast reliable interactive data analysis *//
+//* plot.cpp: plot on-line files *//
+//* (C) Joachim Wuttke 1990-, v2(C++) 2001- *//
+//* http://frida.sourceforge.net *//
+//**************************************************************************//
+
+#include <stdlib.h>
+#include <math.h>
+#include <iostream>
+#include <algorithm>
+
+#include "mystd.h"
+#include "olm.h"
+#include <ask_simple_value.h>
+#include "dualplot.h"
+#include "plot.h"
+
+#include <boost/format.hpp>
+using boost::format;
+
+namespace NPlot {
+/* auxiliary routines */
+ double FindAngleCos( double x0,double y0,
+ double x1,double y1,
+ double x2,double y2 );
+ void FindNextPoint( const COlc *fc,
+ double xtm, double ytm, double dx,
+ double* xt, double* yt, uint k, uint j );
+ void InterY( double xtm,double ytm,double xt,double yt,
+ double* xtt,double *ytt,double border);
+};
+
+double NPlot::FindAngleCos(double x1,double y1, //point 1
+ double x2,double y2, //point 2
+ double x3,double y3 )//point 3
+// Li Huang mar10
+{
+ /*
+ (1) ________(2)
+ \)cos
+ \
+ \(3)
+
+ (x1,y1), (x2,y2)
+ (x2,y2) and (x3,y3)
+ a = x2 - x1, b = y2 - y1
+ c = x3 - x2, d = y3 - y2.
+ The so called dot product yields the equation
+ Cosine( Theta ) = ( (a x c) + (b x d) )/(Sqrt(a^2 + b^2) x Sqrt(c^2 + d^2) ) = Z.
+ */
+ double a,b,c,d,cos;
+ a = x2 - x1;
+ b = y2 - y1;
+ c = x3 - x2;
+ d = y3 - y2;
+ cos = ( (a * c) + (b * d) )/(sqrt(a*a + b*b) * sqrt(c*c + d*d) );
+ return cos;
+}
+
+//! The next point to be plotted is returned as *xt, *yt.
+
+void NPlot::FindNextPoint( const COlc *fc,
+ double xtm, double ytm, double dx,
+ double* xt, double* yt, uint k, uint j )
+// Li Huang mar10
+{
+ uint i = 0;
+ double cos123 = 0;
+ double cos234 = 0;
+ double x1 = xtm;
+ double y1 = ytm;
+ double x2 = *xt;
+ double y2 = *yt;
+ double x3,y3;
+ double x4,y4;
+ /*
+ /)o1__ (3)
+ (2)_____________
+ / \)o2
+ / \
+ (1)/ \(4)
+
+ (xt,yt)__________(xt+dx,yt+dx)
+ / \
+ / \
+ (xtm,ytm)/ \(xt+2*dx,yt+2*dx)
+
+ o1 < 1/4 degree && o2 < 1/4 degree
+ cos(1/4 degree) = 0.99999
+ */
+ do{
+ x3 = x2 + dx;
+ y3 = fc->T->tree_curve_val_sca( x3, k, j );
+ x4 = x3 + dx;
+ y4 = fc->T->tree_curve_val_sca( x4, k, j );
+ cos123 = FindAngleCos(x1,y1,x2,y2,x3,y3);
+ cos234 = FindAngleCos(x2,y2,x3,y3,x4,y4);
+ dx = dx/2;
+ i++;
+ }while ( ( cos123 < 0.99999 || cos234 < 0.99999 ) && i <= 7);
+ *xt = x3;
+ *yt = y3;
+}
+
+void NPlot::InterY(double xtm,double ytm,double xt,double yt,
+ double* xtt,double *ytt,double border)
+{
+ *ytt = border;
+ *xtt = (*ytt-ytm)*(xt-xtm)/(yt-ytm)+xtm;
+}
+
+
+//! Plot spectrum.
+
+void NPlot::Plot( class CPlot *plot, bool add )
+{
+ static CList JSel;
+ CEle *e;
+ uint k, j;
+ const COld *fd;
+ const COlc *fc;
+ CSpec S;
+ vector<double> novec;
+ static int pstyle = 0, cstyle = 0;
+ static string jSel = "";
+ NOlm::JSelAsk( "Plot which spectra", &jSel, &JSel );
+
+ if (!add) {
+
+ pstyle = 0;
+ cstyle = 0;
+
+ // determine x-range ?
+
+ double inf, sup;
+ if (!plot->X.finite()) {
+ NOlm::IterateD fiter;
+ const vector<double> *vx;
+ bool first = true;
+ while((fd=fiter())) {
+ JSel.evaluate( 0, fd->nSpec()-1 );
+ 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 )
+ continue;
+ 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( first )
+ throw string( "x range is empty" );
+ plot->X.setRoundedLimits( "x", inf, sup );
+ }
+
+ // determine y-range ?
+ if (!plot->Y.finite()) {
+ bool first=true;
+ NOlm::IterateEle eiter;
+ int npts = 100;
+ const CSpec *s;
+ while( (e=eiter()) ){
+ k = eiter.SelNo();
+ fc = e->C();
+ if( fc && fc->kconv!=-1 )
+ continue;
+ fd = e->D();
+ JSel.evaluate( 0, e->nSpec()-1 );
+ for( uint iv=0; iv<JSel.size(); ++iv ){
+ uint j = JSel.V[iv];
+ if( fd ){
+ s = &(fd->VS[j]);
+ } else {
+ S.x.resize( npts );
+ for ( uint i=0; i<npts; i++ )
+ S.x[i] = plot->X.plotcoord2value( i/(npts-1.0) );
+ fc->T->tree_curve_val_vec( &(S.y), S.x, k, j );
+ s = &S;
+ }
+ for (uint i=0; i<s->size(); i++) {
+ if( !plot->X.contains( s->x[i] ) )
+ continue;
+ if( plot->Y.logflag && s->y[i]<=0 )
+ continue;
+ if (first) {
+ inf=s->y[i], sup=s->y[i];
+ first = false;
+ } else {
+ if (s->y[i]<inf) inf=s->y[i];
+ if (s->y[i]>sup) sup=s->y[i];
+ }
+ }
+ }
+ }
+ if( first )
+ throw string( "y range is empty" );
+ plot->Y.setRoundedLimits( "y", inf, sup );
+ }
+
+ // build label: (zur Zeit nur vom ersten File; definiere +=)
+ CCoord xCo = NOlm::MOM[*(NOlm::FSel.V.begin())].P()->xco;
+ CCoord yCo = NOlm::MOM[*(NOlm::FSel.V.begin())].P()->yco;
+
+ // draw new frame:
+ plot->clearFrame();
+ plot->plotFrame( xCo.ps_str(), yCo.ps_str() );
+ }
+
+ // plot data:
+ NOlm::IterateEle eiter;
+ while( (e=eiter()) ){
+ k = eiter.SelNo();
+ fc = e->C();
+ fd = e->D();
+ JSel.evaluate( 0, e->nSpec()-1 );
+ for ( uint iv=0; iv<JSel.size(); ++iv){
+ j = JSel.V[iv];
+ // cout << "DEBUG edif plot k="<<k<<", j="<<j<<"\n";
+ if ( fd ) {
+ // plot data
+ const CSpec *s;
+ s = &( fd->VS[j] );
+ uint n=s->x.size();
+ if ( n!=s->y.size() )
+ throw string( "BUG: plot: x.size<>y.size" );
+ // Filter: x,y -> xp,yp if inside frame
+ uint np=0, nxl=0, nxh=0, nyl=0, nyh=0;
+ vector<double> xp, yp, dyp;
+ for ( uint i=0; i<n; i++ ) {
+ double x = s->x[i];
+ if( !plot->X.contains(x) ) {
+ if ( x<=plot->X.inf ){
+ x = plot->X.inf;
+ nxl++;
+ }
+ if ( x>=plot->X.sup ){
+ x = plot->X.sup;
+ nxh++;
+ }
+ if ( !plot->X.force )
+ continue;
+ }
+ double y = s->y[i];
+ if( !plot->Y.contains(y) ) {
+ if ( y<=plot->Y.inf ){
+ y = plot->Y.inf;
+ nyl++;
+ }
+ if ( y>=plot->Y.sup ){
+ y = plot->Y.sup;
+ nyh++;
+ }
+ if ( !plot->Y.force )
+ continue;
+ }
+ if( np==plot->maxpoints && np<n ) {
+ printf("reached maxpoints at %g\n", s->x[i]);
+ }
+ xp.push_back( x );
+ yp.push_back( y );
+ if( s->dy.size() )
+ dyp.push_back( s->dy[i] );
+ np = xp.size();
+ if( np > plot->maxpoints ) // continue in large steps
+ i += plot->maxpoints;
+ }
+ if ( np==0 ) {
+ printf( "file %u spec %u:"
+ " all %u points outside plot range:\n",
+ k, j, n );
+ if( nxl )
+ printf( " %u points < xmin\n", nxl );
+ if( nxh )
+ printf( " %u points > xmax\n", nxh );
+ if( nyl )
+ printf( " %u points < ymin\n", nyl );
+ if( nyh )
+ printf( " %u points > ymax\n", nyh );
+ continue;
+ } else {
+ plot->addSpec( false, pstyle++,
+ xp, yp, dyp, *(e->Z(j)),
+ e->P()->xco.str(), e->P()->yco.str(),
+ "data file " + strg(k) +
+ " spectrum "+strg(j) );
+ }
+ } else if ( fc && fc->kconv!=-1 ) {
+ // plot curve, but only at x points of convolutand
+ uint kconv = fc->kconv;
+ uint jconv = CTree::setConvJ( k, j, kconv );
+ COld *fconv = NOlm::MOM[kconv].D();
+ if ( !fconv )
+ throw string( "convolutand is not a data file" );
+ vector<double> xp, yp;
+ for( uint i=0; i<fconv->VS[jconv].size(); ++i ){
+ double x = fconv->VS[jconv].x[i];
+ if( plot->X.contains( x ) )
+ xp.push_back( x );
+ }
+ fc->T->tree_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] ) &&
+ plot->Y.contains( yp[i] ) ) )
+ continue;
+ xo.push_back( xp[i] );
+ yo.push_back( yp[i] );
+ }
+ if( xo.size()==0 ){
+ cout << "curve k="<<strg(k)<<", j="<<strg(j)<<" is empty\n";
+ return;
+ }
+ plot->addSpec( true, cstyle++, xo, yo, novec, *(e->Z(j)),
+ e->P()->xco.str(), e->P()->yco.str(),
+ "curve file "+strg(k)+" spectrum "+strg(j) );
+ } else if ( fc && false ) {
+ // PRESENTLY UNUSED
+ // plot ordinary curve
+ // straight vectorial version, for rapid plotting
+ int npts = 297;
+ vector<double> xc(npts), yc, xp, yp;
+ for( uint i=0; i<npts; ++i )
+ xc[i] = plot->X.plotcoord2value( i/(npts-1.0) );
+ fc->T->tree_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]) ) {
+ xp.push_back( xc[i] );
+ yp.push_back( yc[i] );
+ }
+ }
+ plot->addSpec( true, cstyle++, xp, yp, novec, *(e->Z(j)),
+ e->P()->xco.str(), e->P()->yco.str(),
+ "curve file "+strg(k)+" spectrum "+strg(j) );
+ } else if ( fc ) {
+ // plot ordinary curve, taking care:
+ // - stop/start line when leaving/reentering the frame (done)
+ // - interpolate when crossing the horizontal border (done)
+ // - interpolate at discontinuities (done)
+ //iamhere
+ uint npts = 97;
+ double dx = plot->X.plotcoord2value( 1/(npts-1.0) ) -
+ plot->X.plotcoord2value( 0/(npts-1.0) );
+ vector< vector<double> > xpv, ypv;
+ vector<double> xp, yp;
+ 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)
+ while ( (!plotend) ) {
+ if( xpv.size()>180 )
+ throw string( "curve oscillates too fast for "
+ "interpolating plot mode" );
+ if ( !plotlast ) {
+ if ( plotfirst ){
+ xt = plot->X.plotcoord2value( 0.0 );
+ yt = fc->T->tree_curve_val_sca( xt, k, j );
+ } else {
+ FindNextPoint( fc, xtmo, ytmo, dx,
+ &xt, &yt, k, j );
+ }
+ }
+ if ( plot->X.contains(xt) ){// if xt <= X.sup
+ if ( plot->Y.contains(yt) ) { //if in
+ if ( plotfirst ) { //a new curve starts
+ //with the first point as the starting point
+ xp.push_back( xt );
+ yp.push_back( yt );
+ plotfirst = false;
+ started = true;
+ }
+ else if ( !started ) { //a new curve starts
+ //i.e. the case: out -> in
+ double xtt,ytt;
+ double border = ytm > plot->Y.sup ?
+ plot->Y.sup : plot->Y.inf;
+ //using the previous point and current point to
+ //interpolate the intersection point
+ InterY(xtm,ytm,xt,yt,&xtt,&ytt,border);
+ xp.push_back( xtt );
+ yp.push_back( ytt );
+ xp.push_back( xt );
+ yp.push_back( yt );
+ started = true;
+ }
+ else if ( started ) { //continues a started curve
+ //i.e. the case in -> in
+ xp.push_back( xt );
+ yp.push_back( yt );
+ }
+ }
+ else if ( started ) {
+ // a curve ends
+ //i.e. the case: in -> in -> out
+ double xtt,ytt;
+ double border = yt > plot->Y.sup ?
+ plot->Y.sup : plot->Y.inf;
+ //using the previous point and current point to
+ //interpolate the intersection point
+ InterY(xtm,ytm,xt,yt,&xtt,&ytt,border);
+ xp.push_back( xtt );
+ yp.push_back( ytt );
+ started = false;
+ ended = true;
+ }
+ else if ( plotfirst ) { //a new curve starts
+ //with the first point is not inside the plotting
+ plotfirst = false;
+ }
+ if ( xt == plot->X.sup ) {
+ ended = true;
+ plotlast = true;
+ plotend = true;
+ }
+ } else if ( ( xt > plot->X.sup ) && !plotlast ){
+ xt = plot->X.sup;
+ plotlast = true;
+ } else
+ throw string( "PROGRAM ERROR plot in fc: no ending" );
+
+ xtmo = xtm;
+ ytmo = ytm;
+ xtm = xt;
+ ytm = yt;
+ if ( ended ){
+ xpv.push_back(xp);
+ ypv.push_back(yp);
+ xp.clear();
+ yp.clear();
+ ended = false;
+ }
+ } // end loop over points
+
+ for ( uint ipv=0; ipv<xpv.size();++ipv) {
+ if(!xpv[ipv].empty()){
+ plot->addSpec( true, cstyle, xpv[ipv], ypv[ipv],
+ novec, *(e->Z(j)),
+ e->P()->xco.str(), e->P()->yco.str(),
+ "curve file "+strg(k)+
+ " spectrum "+strg(j) );
+ }
+ }
+ cstyle++;
+ } else
+ throw string( "PROGRAM ERROR plot invalid *e" );
+ }
+ plot->plotDoclines( e->P()->lDoc );
+ if( fc )
+ plot->plotDoclines( fc->pInfo() );
+ }
+ plot->showSpecs();
+}
diff --git a/pub/src/plot.h b/pub/src/plot.h
new file mode 100644
index 00000000..b6d9e924
--- /dev/null
+++ b/pub/src/plot.h
@@ -0,0 +1,3 @@
+namespace NPlot {
+ void Plot( class CPlot *plot, bool add );
+};
diff --git a/pub/src/plotaux.cpp b/pub/src/plotaux.cpp
deleted file mode 100644
index 8e2e8b4a..00000000
--- a/pub/src/plotaux.cpp
+++ /dev/null
@@ -1,138 +0,0 @@
-//**************************************************************************//
-//* FRIDA: fast reliable interactive data analysis *//
-//* plotaux.cpp: ticks and tacks *//
-//* (C) Joachim Wuttke 1990-, v2(C++) 2001- *//
-//* http://frida.sourceforge.net *//
-//**************************************************************************//
-
-#include <string>
-#include <stdlib.h>
-#include <stdio.h> // for debug
-#include <math.h>
-using namespace std;
-
-#include "plotaux.h"
-#include "mystd.h"
-
-void plotaux::calc_ticks(int logscale, double vmin, double vmax,
- int *ntack, double **tack, int *ntpt, double *ticklim)
-{
- int ir, nt, i;
- double r, rd, rdr, d, vsub, vsup, tack0, tackn, dtack;
-
- if (vmin >= vmax)
- throw "PROGRAM ERROR/ Ticks/ Bad Limits: " +
- strg(vmin) + ", " + strg(vmax);
-
- if (!logscale) {
- r = log10(vmax - vmin);
- ir = (int) ((r<0) ? r-1 : r);
- rd = r - ir; // fractional part of r
- rdr = pow(10., rd);
- if (rdr > 9.99999) {
- dtack = 2.5 * pow(10., ir); // spacing between tacks
- *ntpt = 5; // ticks per tack
- } else if (rdr > 5.00001) {
- dtack = 2 * pow(10., ir);
- *ntpt = 4;
- } else if (rdr > 2.500005) {
- dtack = 1 * pow(10., ir);
- *ntpt = 5;
- } else if (rdr > 1.600003) {
- dtack = .5 * pow(10., ir);
- *ntpt = 5;
- } else if (rdr > 1.250002) {
- dtack = .4 * pow(10., ir);
- *ntpt = 4;
- } else {
- dtack = .25 * pow(10., ir);
- *ntpt = 5;
- }
- d = vmin / dtack;
- tack0 = dtack * (int) (d<=0?d-1:d);
- d = vmax / dtack;
- tackn = dtack * ((int) (d<0?d-1:d) + 1);
- vsub = vmin - (vmax-vmin)*1.e-4;
- vsup = vmax + (vmax-vmin)*1.e-4;
- nt = (int) ((tackn-tack0) / dtack) + 3; // allocate enough
- *tack = (double*) malloc(sizeof(double)*nt);
- *ntack = 0;
- for (i=0; i<nt; i++) {
- d = tack0 + i*dtack;
- if (vsub<=d && d<=vsup)
- (*tack)[(*ntack)++] = d;
- }
-
- ticklim[0] = tack0;
- ticklim[1] = (*tack)[(*ntack)-1] + dtack;
-
- } else { // log scale
-
- static double eins = 1.0000009999999999;
- double rlgmin, rlgmax, rlgrel, d1, d2, r0;
- int minexp, maxexp, incr;
-
- if (vmin <= 0)
- throw string( "PROGRAM ERROR/ Ticks/ negative log argument" );
-
- rlgmin = log10(vmin);
- rlgmax = log10(vmax);
- rlgrel = rlgmax - rlgmin;
- if (rlgrel > 40)
- throw string( "Excessive log range" );
-
- minexp = (rlgmin>0 ?
- (int) (rlgmin+1e-6) :
- ((int) (rlgmin-1e-6) -1));
- maxexp = (rlgmax>0 ?
- (int) (rlgmax+1e-6) :
- ((int) (rlgmax-1e-6) -1));
-
- incr = (int) (rlgrel / 7 + 1); // #decades per tack
- if (incr==2) incr = 1; // preferences: 10^1, 10^3, 10^6 etc.
- if (incr>3) incr = 3 * (incr % 3);
-
- if (incr > 1) {
- minexp -= minexp % incr;
- maxexp -= maxexp % incr;
- }
-
- /* - Set large ticks array : */
- static int mt=30;
- *tack = (double*) malloc(sizeof(double)*mt);
-
- *ntack = 0;
- d1 = vmin / eins;
- d2 = vmax * eins;
- for ( i = minexp; incr < 0 ? i >= maxexp : i <= maxexp; i += incr ) {
- r0 = pow(10., i);
- if (d1<=r0 && r0<=d2) {
- if(*ntack>=mt)
- throw string( "PROG ERR too many tacks" );
- (*tack)[(*ntack)++] = r0;
- }
- }
-
- if (*ntack >= 1) {
- ticklim[0] = (*tack)[0] / pow(10., incr);
- ticklim[1] = (*tack)[(*ntack)-1] * pow(10., incr);
- } else {
- // untested. just to allow plots when no tacks in window
- ticklim[0] = pow( 10., minexp );
- ticklim[1] = pow( 10., maxexp+1 );
- }
-
- /** determine # ticks / tack **/
- if (incr==1 && rlgrel<=7) {
- *ntpt = 9;
- } else if (incr==1 && rlgrel<=7) {
- *ntpt = 3; // 1-2-5-10 - Schritte
- } else if (incr==1) {
- *ntpt = 1;
- } else {
- *ntpt = -incr;
- // means incr ticks per tack
- // means 1 tick per decade
- }
- }
-}
diff --git a/pub/src/plotaux.h b/pub/src/plotaux.h
deleted file mode 100644
index a66a336e..00000000
--- a/pub/src/plotaux.h
+++ /dev/null
@@ -1,6 +0,0 @@
-// This is plotaux.h. For details, see plotaux.cpp
-
-namespace plotaux {
- void calc_ticks(int logscale, double vmin, double vmax,
- int *ntack, double **tack, int *ntpt, double *ticklim);
-}
diff --git a/pub/src/rssm.cpp b/pub/src/rssm.cpp
index a714101a..376080de 100644
--- a/pub/src/rssm.cpp
+++ b/pub/src/rssm.cpp
@@ -9,6 +9,7 @@
#include <stdio.h>
#include <time.h>
#include <iostream>
+#include <fstream>
#include "mystd.h"
#include "olm.h"
--
GitLab