dualplot.cpp

//**************************************************************************//
//* FRIDA: flexible rapid interactive data analysis                        *//
//* dualplot.cpp: different mechanisms for screen and paper output         *//
//* (C) Joachim Wuttke 1990-, v2(C++) 2001-                                *//
//* http://www.messen-und-deuten.de/frida                                  *//
//**************************************************************************//

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <iostream>
#include <math.h>
#include <fcntl.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <boost/format.hpp>

#include "mystd.h"
#include "dualplot.h"

using boost::format;


//! Constructor for plot window: setup for gnuplot and postscript.

CPlot::CPlot( uint _iPlot, bool _logx, bool _logy ) :
    iPlot( _iPlot ), X( _logx ), Y( _logy ),
    maxpoints(20000), with_errors(true), equipoints(7), refine(true)
{
    // == Initialization for Gnuplot ==

    // Create a named pipe (FIFO) that will transmit our commands to Gnuplot.
    string fn_fifo = string("/tmp/gnuplot-") + getenv( "LOGNAME" );
    mystd::system( "rm -f "+fn_fifo );
    if (mkfifo( fn_fifo.c_str(), 0666 ))
        throw "SYSTEM ERROR cannot make fifo " + fn_fifo +
            ": will not be able to print";

    // Start a Gnuplot that will read from our pipe.
    mystd::system( "gnuplot -title " + strg(iPlot) + " -noraise" +
                   " < " + fn_fifo + " &" );

    // Open our pipe so that we can write commands to it
    //   (we use 'open' instead of 'fopen' or 'ofstream',
    //    because we need non-blocking mode).
    if ( ! ( gp_fifo = open( fn_fifo.c_str(), O_WRONLY ) ) )
        throw "SYSTEM ERROR cannot open fifo " + fn_fifo +
            ": will not be able to print";
    fcntl( gp_fifo, F_SETFL, O_NONBLOCK );

    // Now the initialization _within_ Gnuplot.
    gp_write( string("set terminal x11") );

    // == Initialization for PostScript ==
    ps_fnum=0;
}


//! Clear plot frame and buffer.

void CPlot::clearFrame()
{
    // clear gnuplot tmp files
    gp_fno = 0;
    gp_fnames = "";

    // reset buffers for postscript output
    string cmd;
    ps_accu.clear();
    ps_accu.push_back( "\n%% output created by frida2\n");
    ps_snum = 0;
    ps_pnum = 0;
    ps_cnum = 0;
    ps_Doc.clear();
}


//! Change one setup parameter per command.

void CPlot::setAux( string cmd )
{
    if ( cmd=="gxl" ) {
        X.setLog( !X.logflag );
        printf( "set x %s\n", X.logflag ? "log" : "lin" ); }
    else if ( cmd=="gxl+" )
        X.setLog( true );
    else if ( cmd=="gxl-" )
        X.setLog( false );
    else if ( cmd=="gyl" ) {
        Y.setLog( !Y.logflag );
        printf( "set y %s\n", Y.logflag ? "log" : "lin" ); }
    else if ( cmd=="gyl+" )
        Y.setLog( true );
    else if ( cmd=="gyl-" )
        Y.setLog( false );
    else if ( cmd=="gxf" ) {
        X.force = !X.force;
        printf( "force x %s\n", X.force ? "on" : "off" ); }
    else if ( cmd=="gxf+" )
        X.force = true;
    else if ( cmd=="gxf-" )
        X.force = false;
    else if ( cmd=="gyf" ) {
        Y.force = !Y.force;
        printf( "force y %s\n", Y.force ? "on" : "off" ); }
    else if ( cmd=="gyf+" )
        Y.force = true;
    else if ( cmd=="gyf-" )
        Y.force = false;
    else if ( cmd=="ge" )
        with_errors = !with_errors;
    else if ( cmd=="ge+" )
        with_errors = true;
    else if ( cmd=="ge-" )
        with_errors = false;
    else
        throw "unknown command " + cmd;
}


//! Plot coordinate frame (axes, ticks, labels).

void CPlot::plotFrame( string xlabel, string ylabel )
{
    gp_write( "set nologscale" );
    string whichlog="";
    if (X.logflag) whichlog += "x";
    if (Y.logflag) whichlog += "y";
    if (whichlog!="")
        gp_write( "set logscale " + whichlog );

    // wups:
    snprintf( outlin, mLin, "\n%d %g %g xSetCoord\n", 
              X.logflag, X.inf, X.sup );
    ps_accu.push_back( outlin );
    snprintf( outlin, mLin, "%d %g %g ySetCoord\n", 
              Y.logflag, Y.inf, Y.sup );
    ps_accu.push_back( outlin );
    snprintf( outlin, mLin, "%% %d %g %g %d zSetCoord\n\n",
              0, 0., 0., 0 );
    ps_accu.push_back( outlin );

    int ntack, ntpt;
    double *tack, ticklim[2];

    ps_accu.push_back( "\n/xPlotFrame {\n" );
    if ( X.logflag && X.inf<= 0 )
        throw "BUG: x log incompatible with limits " + X.str();
    X.calc_ticks( &ntack, &tack, &ntpt, ticklim );
    ps_ticktack(ntack, tack, ntpt, ticklim, &X);
    free(tack);
    snprintf( outlin, mLin-4, "  {(%s", xlabel.c_str() );
    strncat( outlin, ")}\n", mLin );
    ps_accu.push_back( outlin );
    ps_accu.push_back( "   0 10   0  0     0  90 "
                       "OneAxx Axx Tic Tac xNumL %% low x axis\n" );
    ps_accu.push_back( "   0 10   0 10     0 270 "
                       "OneAxx Axx Tic Tac       %% top x axis\n" );
    ps_accu.push_back( "  xCL\n" );
    ps_accu.push_back( "} def\n" );
    
    ps_accu.push_back( "\n/yPlotFrame {\n" );
    if ( Y.logflag && Y.inf<= 0 )
        throw "BUG: y log incompatible with limits " + Y.str();
    Y.calc_ticks( &ntack, &tack, &ntpt, ticklim );
    ps_ticktack(ntack, tack, ntpt, ticklim, &Y);
    free(tack);
    snprintf( outlin, mLin, "   {(%s", ylabel.c_str() );
    strncat( outlin, ")}\n", mLin );
    ps_accu.push_back( outlin );
    ps_accu.push_back( "   0 10   0  0    90   0 "
                       "OneAxx Axx Tic Tac yNumL %% left y axis\n" );
    ps_accu.push_back( "   0 10  10  0    90 180 "
                       "OneAxx Axx Tic Tac       %% right y axis\n" );
    ps_accu.push_back( "  yCL\n" );
    ps_accu.push_back( "} def\n" );
    ps_accu.push_back( "\n%% modeDD\nplotbefore\n" );
}


//! Plot one spectrum.

void CPlot::addSpec( bool as_line, int style_no,
                     const vector<double>& xp,
                     const vector<double>& yp, const vector<double>& dyp,
                     const vector<double>& z,
                     string xco, string yco, string info
    )
{
    static const int mColor = 6;
    static int color[mColor] = { 0x880000, 0x008800, 0x000088,
                                 0x006666, 0x660066, 0x666600 };
    // Checks:
    uint np=xp.size();
    if ( !np )
        throw string( "BUG: NPLot::Line x.size=0" );
    if ( np!=yp.size() )
        throw string( "BUG: NPLot::Line x.size<>y.size" );

    // Prepare for live display, to be shown by showSpecs():
    string gp_fnam = str( format( "/tmp/%s-%d-%03d.gnu" )
                          % getenv("LOGNAME") % iPlot % gp_fno++ );
    if (gp_fnames!="") gp_fnames += ", ";
    gp_fnames += string("\"") + gp_fnam + "\" notitle";
    if ( as_line )
        gp_fnames += str( format( " with lines lt 1 lc rgb \"#%6x\"" )
                          % color[ style_no % mColor ] );
    else if( with_errors && dyp.size() )
        gp_fnames += " with errorbars";
    FILE *gp_fd;
    if (!(gp_fd = fopen(gp_fnam.c_str(), "w")))
        throw "cannot save gnuplot data to " + gp_fnam;
    uint nout = 0;
    for (uint i=0; i<np; i++){
        if( isinf(xp[i]) || isinf(yp[i]) )
            throw "Data point number " + strg(i) + " is invalid: x=" +
                strg(xp[i]) + ", y=" + strg(yp[i]);
        if( xp[i]<X.inf || xp[i]>X.sup )
            throw "Plot::Line: x["+strg(i)+"]="+strg(xp[i])+" out of range";
        if( yp[i]<Y.inf || yp[i]>Y.sup )
            throw "Plot::Line: y["+strg(i)+"]="+strg(yp[i])+" out of range";
        if( with_errors && dyp.size() )
            fprintf(gp_fd, "%16.8g %16.8g %16.8g\n", xp[i], yp[i], dyp[i] );
        else
            fprintf(gp_fd, "%16.8g %16.8g\n", xp[i], yp[i]);
        nout++;
    }
    fclose(gp_fd);
    if( !nout )
        throw "no points in frame: " + info;


    // Postscript copy:
    snprintf( outlin, mLin, "\n%3u [", ++ps_snum );
    ps_accu.push_back( outlin );
    for (uint i=0; i<z.size(); i++){
        snprintf( outlin, mLin, " %12g", z[i]);
        ps_accu.push_back( outlin );
    }
    snprintf( outlin, mLin, " ] zValues\n" );
    ps_accu.push_back( outlin );
    if ( as_line )
        snprintf( outlin, mLin, "%2d cstyle", 1 );
    else
        snprintf( outlin, mLin, "%2d pstyle", style_no+1 );
    ps_accu.push_back( outlin );
    snprintf( outlin, mLin-2, " %% (%s -> %s)",
              xco.c_str(), yco.c_str() );
    strncat( outlin, "\n", mLin );
    ps_accu.push_back( outlin );
    for (uint i=0; i<np; i++) {
        snprintf( outlin, mLin,
                  "%6.3f %6.3f %6.3f t%c %% %13.7g wx %13.7g wy\n",
                  X.pc(xp[i]), Y.pc(yp[i]), 
                  dyp.size() ? Y.pcerr(yp[i],dyp[i]) : 0,
                  i==0 ? 'i' : i==np-1 ? 'f' : ' ',
                  xp[i], yp[i] );
        ps_accu.push_back( outlin );
    }
}

//! Live display as prepared by addSpec.

void CPlot::showSpecs()
{
    if ( gp_fnames!="" )
        gp_write( "plot " + 
                  str( format( "[%12.8g:%12.8g] [%12.8g:%12.8g] " )
                       % X.inf % X.sup % Y.inf % Y.sup ) +
                  gp_fnames );
}


//! Add documentation lines to postscript output.

void CPlot::plotDoclines( const vector<string>& lDoc )
{
    for (uint i=0; i<lDoc.size(); i++)
        ps_Doc.push_back(lDoc[i]);
}


//! Write buffered plot to postscript file.

void CPlot::writePostscript( string ps_outdir, string ps_head, string ps_dict )
{
    // construct output file name:
    FILE *pssav;
    string cmd, outf;
    while(1) {
        if (ps_fnum>=999)
            throw "graph file number overflow";
        outf = mystd::wordexp_unique( 
            ps_outdir + str( format( "l%d" ) % ++ps_fnum ) + "." +
            ( ps_dict=="" ? "psa" : "ps" ) );
        if( !mystd::file_exists( outf.c_str() ) )
            break; // legal exit
    }
    cout << "save plot in " << outf << "\n";

    // copy headers to output file:
    cmd = string("cat ") +
        ps_dict + " " + // ps_dict may be ""
        ps_head + " > " +
        outf + "\n";
    mystd::system( cmd );
    
    // append specific output to output file:
    if ( !(pssav = fopen( outf.c_str(), "a+" )) )
        throw "cannot append contents to file " + outf;
    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 );
    }

    // additional output (do not append this to ps_accu to allow
    // further incrementation of ps_accu):
    fprintf( pssav, "\nblack 0 -4 13 1.65 NewList\n" );
    for ( uint i=0; i<ps_Doc.size(); i++ ) 
        fprintf( pssav, "{(%s)} infline\n", ps_Doc[i].c_str() );
	
    fprintf( pssav, 
             "\n{(%s)}  /filename exch def 10 -2.8 18 showfilename\n\n"
            " EndFrame\n", outf.c_str() );

    // output completed:
    fclose(pssav);
}


//! Info line to characterize this plot window.

string CPlot::info() const
{
    string ret;
    ret  =   "x: " + X.info();
    ret += "  y: " + Y.info();
    return ret;
}


//! Send one line to gnuplot fifo.

void CPlot::gp_write( string in )
{
    string out = in + "\n";
    // cout << "monitor gnuplot driver: '" << out << "'\n";
    if( write( gp_fifo, out.c_str(), out.size() ) <= 0 )
        throw string( "could not write to gp_fifo" );
}


//! Format ticks and tacks for postscript file.

void CPlot::ps_ticktack(
    uint ntack, double *tack, int ntpt, double *ticklim, CAxis *A )
{
    uint i;
    ps_accu.push_back( "  [\n" );
    if (A->logflag && ( tack[0]<1e-3 || tack[ntack-1]>1e3 )) {
        for (i=0; i<ntack; i++) {
            snprintf( outlin, mLin,
                      "   %9.6f {(10)(%d)sp()} %%{(%g)}\n", 
                      A->pc(tack[i]), (int)(log10(tack[i])),
                      (float) tack[i]);
            ps_accu.push_back( outlin );
        }
    } else {
        for (i=0; i<ntack; i++) {
            snprintf( outlin, mLin, "   %9.6f {(%g)}\n", 
                      A->pc(tack[i]), (float) tack[i]);
            ps_accu.push_back( outlin );
        }
    }
    ps_accu.push_back( "  ] SetTacVec\n" );
    snprintf( outlin, mLin, "  %g %g %d %d SetTicVec%s\n", 
              A->pc(ticklim[0]), A->pc(ticklim[1]), ntack+2, ntpt,
              (A->logflag? "Log" : "Lin"));
    ps_accu.push_back( outlin );
}



namespace NPloWin {
    vector<CPlot*> Plots;
    uint nPlot, iPlot;

    void initialize();
};

//! Initialize plot windows.

void NPloWin::initialize()
{
    using namespace NPloWin;
    nPlot = 0; // shorthand for Plots.size()
    Plots.push_back( new class CPlot( nPlot++, false, false ) );
    Plots.push_back( new class CPlot( nPlot++, true,  false ) );
    Plots.push_back( new class CPlot( nPlot++, false, true ) );
    Plots.push_back( new class CPlot( nPlot++, true,  true ) );
    iPlot = 0; // current plot window
}