binaries/src/disembl/Tisean_3.0.1/source_c/poincare.c

   1 /*
   2  *   This file is part of TISEAN
   3  *
   4  *   Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
   5  *
   6  *   TISEAN is free software; you can redistribute it and/or modify
   7  *   it under the terms of the GNU General Public License as published by
   8  *   the Free Software Foundation; either version 2 of the License, or
   9  *   (at your option) any later version.
  10  *
  11  *   TISEAN is distributed in the hope that it will be useful,
  12  *   but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  *   GNU General Public License for more details.
  15  *
  16  *   You should have received a copy of the GNU General Public License
  17  *   along with TISEAN; if not, write to the Free Software
  18  *   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  19  */
  20 /*Author: Rainer Hegger, Last modified: Mar 20, 1999 */
  21 #include <string.h>
  22 #include <stdio.h>
  23 #include <stdlib.h>
  24 #include <math.h>
  25 #include <limits.h>
  26 #include "routines/tsa.h"
  27
  28 #define WID_STR "Make a Poincare section"
  29
  30 char *outfile=NULL,dimset=0,compset=0,whereset=0,stdo=1;
  31 char *infile=NULL;
  32 unsigned long length=ULONG_MAX,count,exclude=0;
  33 int dim=2,comp=2,delay=1,dir=0;
  34 unsigned int column=1;
  35 unsigned int verbosity=0xff;
  36 double *series,min,max,average=0.0,where;
  37
  38 void show_options(char *progname)
  39 {
  40   what_i_do(progname,WID_STR);
  41   fprintf(stderr," Usage: %s [Options]\n\n",progname);
  42   fprintf(stderr," Options:\n");
  43   fprintf(stderr,"Everything not being a valid option will be interpreted"
  44           " as a possible"
  45           " datafile.\nIf no datafile is given stdin is read. Just - also"
  46           " means stdin\n");
  47   fprintf(stderr,"\t-l # of points to be used [Default: whole file]\n");
  48   fprintf(stderr,"\t-x #of lines to be ignored [Default: 0]\n");
  49   fprintf(stderr,"\t-c column to read  [Default: 1]\n");
  50   fprintf(stderr,"\t-m embedding dimension [Default: 2]\n");
  51   fprintf(stderr,"\t-d delay [Default: 1]\n");
  52   fprintf(stderr,"\t-q component to cut [Default: last]\n");
  53   fprintf(stderr,"\t-C direction of the cut (0: from below,1: from above)"
  54           "\n\t\t[Default: 0]\n");
  55   fprintf(stderr,"\t-a set crossing at [Default: average of data]\n");
  56   fprintf(stderr,"\t-o outfile [Default: 'datafile'.poin]\n");
  57   fprintf(stderr,"\t-V verbosity level [Default: 1]\n\t\t"
  58           "0='only panic messages'\n\t\t"
  59           "1='+ input/output messages'\n");
  60   fprintf(stderr,"\t-h  show these options\n");
  61   fprintf(stderr,"\n");
  62   exit(0);
  63 }
  64
  65
  66 void scan_options(int n,char** in)
  67 {
  68   char *out;
  69
  70   if ((out=check_option(in,n,'l','u')) != NULL)
  71     sscanf(out,"%lu",&length);
  72   if ((out=check_option(in,n,'x','u')) != NULL)
  73     sscanf(out,"%lu",&exclude);
  74   if ((out=check_option(in,n,'c','u')) != NULL)
  75     sscanf(out,"%u",&column);
  76   if ((out=check_option(in,n,'m','u')) != NULL) {
  77     dimset=1;
  78     sscanf(out,"%u",&dim);
  79   }
  80   if ((out=check_option(in,n,'d','u')) != NULL)
  81     sscanf(out,"%u",&delay);
  82   if ((out=check_option(in,n,'q','u')) != NULL) {
  83     compset=1;
  84     sscanf(out,"%u",&comp);
  85   }
  86   if ((out=check_option(in,n,'C','u')) != NULL)
  87     sscanf(out,"%u",&dir);
  88   if ((out=check_option(in,n,'V','u')) != NULL)
  89     sscanf(out,"%u",&verbosity);
  90   if ((out=check_option(in,n,'a','f')) != NULL) {
  91     whereset=1;
  92     sscanf(out,"%lf",&where);
  93   }
  94   if ((out=check_option(in,n,'o','o')) != NULL) {
  95     stdo=0;
  96     if (strlen(out) > 0)
  97       outfile=out;
  98   }
  99 }
 100
 101 void poincare(void)
 102 {
 103   unsigned long i;
 104   long j,jd;
 105   double delta,xcut;
 106   double time=0.0,lasttime=0.0;
 107   FILE *fout=NULL;
 108
 109   if (!stdo) {
 110     fout=fopen(outfile,"w");
 111     if (verbosity&VER_INPUT)
 112       fprintf(stderr,"Opened %s for writing\n",outfile);
 113   }
 114   else {
 115     if (verbosity&VER_INPUT)
 116       fprintf(stderr,"Writing to stdout\n");
 117   }
 118
 119   if (dir == 0) {
 120     for (i=(comp-1)*delay;i<length-(dim-comp)*delay-1;i++) {
 121       if ((series[i] < where) && (series[i+1] >= where)) {
 122         delta=(series[i]-where)/(series[i]-series[i+1]);
 123         time=(double)i+delta;
 124         if (lasttime > 0.0) {
 125           for (j= -(comp-1);j<=dim-comp;j++) {
 126             if (j != 0) {
 127               jd=i+j*delay;
 128               xcut=series[jd]+delta*(series[jd+1]-series[jd]);
 129               if (!stdo)
 130                 fprintf(fout,"%e ",xcut);
 131               else
 132                 fprintf(stdout,"%e ",xcut);
 133             }
 134           }
 135           if (!stdo) {
 136             fprintf(fout,"%e\n",time-lasttime);
 137             fflush(fout);
 138           }
 139           else {
 140             fprintf(stdout,"%e\n",time-lasttime);
 141             fflush(stdout);
 142           }
 143           count++;
 144         }
 145         lasttime=time;
 146       }
 147     }
 148   }
 149   else {
 150     for (i=(comp-1)*delay;i<length-(dim-comp)*delay-1;i++) {
 151       if ((series[i] > where) && (series[i+1] <= where)) {
 152         delta=(series[i]-where)/(series[i]-series[i+1]);
 153         time=(double)i+delta;
 154         if (lasttime > 0.0) {
 155           for (j= -(comp-1);j<=dim-comp;j++) {
 156             if (j != 0) {
 157               jd=i+j*delay;
 158               xcut=series[jd]+delta*(series[jd+1]-series[jd]);
 159               if (!stdo)
 160                 fprintf(fout,"%e ",xcut);
 161               else
 162                 fprintf(stdout,"%e ",xcut);
 163             }
 164           }
 165           if (!stdo) {
 166             fprintf(fout,"%e\n",time-lasttime);
 167             fflush(fout);
 168           }
 169           else {
 170             fprintf(stdout,"%e\n",time-lasttime);
 171             fflush(stdout);
 172           }
 173           count++;
 174         }
 175         lasttime=time;
 176       }
 177     }
 178   }
 179   if (!stdo)
 180     fclose(fout);
 181 }
 182
 183 int main(int argc,char** argv)
 184 {
 185   char stdi=0;
 186   long i;
 187   double var;
 188
 189   if (scan_help(argc,argv))
 190     show_options(argv[0]);
 191
 192   scan_options(argc,argv);
 193 #ifndef OMIT_WHAT_I_DO
 194   if (verbosity&VER_INPUT)
 195     what_i_do(argv[0],WID_STR);
 196 #endif
 197
 198   infile=search_datafile(argc,argv,&column,verbosity);
 199   if (infile == NULL)
 200     stdi=1;
 201
 202   if (outfile == NULL) {
 203     if (!stdi) {
 204       check_alloc(outfile=(char*)calloc(strlen(infile)+6,(size_t)1));
 205       strcpy(outfile,infile);
 206       strcat(outfile,".poin");
 207     }
 208     else {
 209       check_alloc(outfile=(char*)calloc((size_t)11,(size_t)1));
 210       strcpy(outfile,"stdin.poin");
 211     }
 212   }
 213   if (!stdo)
 214     test_outfile(outfile);
 215
 216   series=(double*)get_series(infile,&length,exclude,column,verbosity);
 217   variance(series,length,&average,&var);
 218   min=max=series[0];
 219   for (i=1;i<length;i++) {
 220     if (series[i] < min) min=series[i];
 221     if (series[i] > max) max=series[i];
 222   }
 223
 224   if (!whereset)
 225     where=average;
 226   if (dimset && !compset)
 227     comp=dim;
 228
 229   if (comp > dim) {
 230     fprintf(stderr,"Component to cut is larger than dimension. Exiting!\n");
 231     exit(POINCARE__WRONG_COMPONENT);
 232   }
 233   if ((where < min) || (where > max)) {
 234     fprintf(stderr,"You want to cut outside the data interval which is [%e,"
 235             "%e]\n",min,max);
 236     exit(POINCARE__OUTSIDE_REGION);
 237   }
 238   poincare();
 239
 240   return 0;
 241 }