-
Notifications
You must be signed in to change notification settings - Fork 0
/
fileplot.cpp
788 lines (755 loc) · 38.5 KB
/
fileplot.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
/*============================================================================
* Daniel J. Greenhoe
*============================================================================*/
/*=====================================
* headers
*=====================================*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include "main.h"
#include "symseq.h"
#include "r1.h"
#include "r2.h"
#include "r3.h"
#include "r4.h"
#include "r6.h"
#include "c1.h"
#include "c6.h"
#include "elliptic.h"
#include "larc.h"
#include "mca.h"
#include "euclid.h"
#include "die.h"
#include "realdie.h"
#include "fairdie.h"
#include "spinner.h"
#include "dna.h"
#include "dnan.h"
#include "fileplot.h"
/*-------------------------------------------------------------------------
* open plot file <filename> and return pointer to file
*-------------------------------------------------------------------------*/
FILE *plot_open(const char *basefilename, const time_t time1, const char *comment){
struct tm *gmt;
gmt = gmtime(&time1);
char tbuffer[80];
char filename[256];
if(strchr(basefilename,'.')==NULL) sprintf(filename,"%s.tex",basefilename);
else sprintf(filename,"%s",basefilename);
FILE *fptr=fopen(filename,"w");
if(fptr==NULL){
fprintf(stderr,"ERROR attempting to open file %s\n",filename);
}
else{
strftime(tbuffer,80,"%Y %B %d %A %I:%M:%S %p UTC",gmt);
fprintf(fptr,"%%========================================================================\n");
fprintf(fptr,"%% Daniel J. Greenhoe \n");
fprintf(fptr,"%% %s\n",tbuffer);
if(strlen(comment)>0)fprintf(fptr,"%% %s\n",comment);
fprintf(fptr,"%% XeLaTeX graphics file for plotting sequence\n");
fprintf(fptr,"%% using the pstricks package(s).\n");
fprintf(fptr,"%% References: http://www.ctan.org/pkg/pstricks-base\n");
fprintf(fptr,"%% http://www.ctan.org/pkg/pstricks-add\n");
fprintf(fptr,"%% http://www.ctan.org/pkg/preview\n");
fprintf(fptr,"%% To produce a pdf file from this file, enter this at the command line:\n");
fprintf(fptr,"%% xelatex <filename>\n");
fprintf(fptr,"%% This file has been generated by the C++ open source program \n");
fprintf(fptr,"%% ssp.exe by Daniel J. Greenhoe \n");
fprintf(fptr,"%%========================================================================\n");
}
return fptr;
}
/*-------------------------------------------------------------------------
* open plot file <filename> and return pointer to file
*-------------------------------------------------------------------------*/
FILE *log_open(const char *basefilename, const time_t time1, const char *comment){
struct tm *gmt;
char tbuffer[80];
char filename[256];
gmt = gmtime(&time1);
strftime(tbuffer,80,"%Y%m%d_%H%M%S",gmt);
if (strstr(basefilename,".log")!=NULL) sprintf(filename,"%s",basefilename);
else if(strstr(basefilename,".xlg")!=NULL) sprintf(filename,"%s",basefilename);
else sprintf(filename,"%s.xlg",basefilename);
printf("open experiment log file \"%s\"...",filename);
FILE *fptr=fopen(filename,"w");
if(fptr==NULL){
fprintf(stderr,"ERROR attempting to open experiment log file %s\n",filename);
exit(EXIT_FAILURE);
}
printf("success.\n");
strftime(tbuffer,80,"%Y %B %d %A %I:%M:%S %p UTC",gmt);
fprintf(fptr,"%%========================================================================\n");
fprintf(fptr,"%% Daniel J. Greenhoe \n");
fprintf(fptr,"%% %s\n",tbuffer);
fprintf(fptr,"%% experiment log file.\n");
if(strlen(comment)>0)fprintf(fptr,"%% %s\n",comment);
fprintf(fptr,"%% This file has been generated by the C++ open source program \n");
fprintf(fptr,"%% ssp.exe by Daniel J. Greenhoe \n");
fprintf(fptr,"%%========================================================================\n");
return fptr;
}
/*-------------------------------------------------------------------------
* open plot file <filename> and return pointer to file
*-------------------------------------------------------------------------*/
int plot_close(FILE *fptr, const time_t time1){
time_t time2;time(&time2);
double seconds;
struct tm *gmt;
gmt = gmtime(&time2);
char tbuffer[80];
strftime(tbuffer,80,"%Y %B %d %A %I:%M:%S %p UTC",gmt);
seconds = difftime(time2,time1);
fprintf(fptr,"%%========================================================================\n");
fprintf(fptr,"%% end \n");
fprintf(fptr,"%% %s\n",tbuffer);
fprintf(fptr,"%% (%.0lf seconds ellapsed)\n",seconds);
fprintf(fptr,"%%========================================================================\n");
fclose(fptr);
return 0;
}
/*-------------------------------------------------------------------------
* print header for XeLaTeX plot files
*-------------------------------------------------------------------------*/
int plot_header(FILE *fptr, const time_t time1, const char *comment){
struct tm *gmt;
gmt = gmtime(&time1);
char tbuffer[80];
strftime(tbuffer,80,"%Y %B %d %A %I:%M:%S %p UTC",gmt);
fprintf(fptr,"%%========================================================================\n");
fprintf(fptr,"%% Daniel J. Greenhoe \n");
fprintf(fptr,"%% %s\n",tbuffer);
if(strlen(comment)>0)fprintf(fptr,"%% %s\n",comment);
fprintf(fptr,"%% XeLaTeX graphics file for plotting sequence\n");
fprintf(fptr,"%% using the pstricks package(s).\n");
fprintf(fptr,"%% References: http://www.ctan.org/pkg/pstricks-base\n");
fprintf(fptr,"%% http://www.ctan.org/pkg/pstricks-add\n");
fprintf(fptr,"%% http://www.ctan.org/pkg/preview\n");
fprintf(fptr,"%% To produce a pdf file from this file, enter this at the command line:\n");
fprintf(fptr,"%% xelatex <filename>\n");
fprintf(fptr,"%% This file has been generated by the C++ open source program \n");
fprintf(fptr,"%% ssp.exe by Daniel J. Greenhoe \n");
fprintf(fptr,"%%========================================================================\n");
return 0;
}
/*-------------------------------------------------------------------------
* print footer for XeLaTeX plot files and log files
*-------------------------------------------------------------------------*/
int plot_footer(FILE *fptr, const time_t time1){
struct tm *gmt;
time_t time2;
double seconds;
char tbuffer[80];
time(&time2);
gmt = gmtime(&time2);
strftime(tbuffer,80,"%Y %B %d %A %I:%M:%S %p UTC",gmt);
seconds = difftime(time2,time1);
fprintf(fptr,"%%========================================================================\n");
fprintf(fptr,"%% end data \n");
fprintf(fptr,"%% %s\n",tbuffer);
fprintf(fptr,"%% (%.0lf seconds ellapsed)\n",seconds);
fprintf(fptr,"%%========================================================================\n");
return 0;
}
/*-------------------------------------------------------------------------
* plot OCS sequence
*-------------------------------------------------------------------------*/
int plot_ocs_seq(const symseq *x, const long start, const long finish, const time_t time1, const char *plottype, const char *filename, const char *comment,FILE *lptr){
time_t time2;
struct tm *gmt;
double seconds;
int i,mode;
long n,m;
char symbol;
FILE *fptr;
char cmdstr[8]="lolli";//"lolli"=default command string
char tbuffer[80];// time buffer
char buf[256];
const char plottypelist[][16]={"lolli","lollipop", "die","dice", "spin","spinner", "dna","genome", "dnan","genomen"};
const int modelist[] ={ 0, 0, 1, 1, 2, 2, 3, 3, 4, 4};
const char cmdstrlist[][16] ={"lolli", "die", "spin", "dna", "dnan"};
const char symlist[][16] ={"0123456", ".ABCDEF", ".ABCDEF", ".ATCG", "NATCG"};
const int Ncmds = sizeof(plottypelist)/sizeof(char)/16;
int Nsym; //number of symbols in a given mode
/*----------------------------------------
* process plottype parameter
*----------------------------------------*/
for(i=0,mode=-1;i<Ncmds;i++)
if(strcmp(plottype,plottypelist[i])==0){
mode=modelist[i];
sprintf(cmdstr,cmdstrlist[mode]);
}
if(mode==-1){
mode=modelist[0];
sprintf(cmdstr,cmdstrlist[mode]);
fprintf(stderr,"\nERROR using plot_ocs_seq(...): unknown plottype string \"%s\".\n",plottype);
exit(EXIT_FAILURE);
}
Nsym = strlen(symlist[mode]);
/*----------------------------------------
* generate TeX plot file
*----------------------------------------*/
if(filename==NULL) fptr=stdout;
else fptr=fopen(filename,"w");
if(fptr==NULL){
sprintf(buf,"\nERROR using plot_die_seq(...): Cannot open output file %s .\n",filename);
printeft(lptr,buf);
exit(EXIT_FAILURE);
}
if(plot_header(fptr, time1, comment))return -1;
fprintf(fptr,"\\input{shelltop.tex}%%\n");
fprintf(fptr,"\\begin{document}%%\n");
fprintf(fptr," \\gsize%%\n");
fprintf(fptr," \\psset{xunit=3mm,yunit=3mm}%%\n");
fprintf(fptr," \\begin{pspicture}(-2,-1.4)(52,%d)%%\n",Nsym);
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," %% axes\n");
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," \\psaxes[linecolor=axis,yAxis=false,showorigin=false,ticks=none,Dx=5,labels=none]{->}(0,0)(0,0)(52,%d)%% x-axis\n",Nsym);
fprintf(fptr," \\psaxes[linecolor=axis,xAxis=false,showorigin=false,labels=none]{->}(0,0)(0,0)(52,%d)%% y-axis\n",Nsym);
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," %% labels along y-axis\n");
fprintf(fptr," %%----------------------------------------\n");
for(i=0;i<Nsym;i++){
symbol = symlist[mode][i];
if(symbol!='.')fprintf(fptr," \\uput[180]{0}(0,%d){$\\%s%c$}%%\n",i,cmdstr,symbol);
}
m = (start==0)? 0 : 2;
if(m)fprintf(fptr," \\rput(1,3.5){$\\cdots$}%%\n");
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," %% labels along x-axis\n");
fprintf(fptr," %%----------------------------------------\n");
for(n=start;n<=finish;n+=5,m+=5)fprintf(fptr," \\uput{2pt}[-90]{0}(%2ld,0){$%ld$}%%\n",m,n);
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," %% data points\n");
fprintf(fptr," %%----------------------------------------\n");
m = (start==0)? 0 : 2;
for(n=start;n<=finish;n++,m++){
symbol=x->get(n);
fprintf(fptr," \\%splot%c{%02ld}%% x[%02ld]=%c\n",cmdstr,symbol,m,n,symbol);
}
fprintf(fptr," \\rput(51,2){$\\cdots$}%%\n");
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," %% end data\n");
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," \\end{pspicture}%%\n");
fprintf(fptr,"\\end{document}%%\n");
time(&time2);
gmt = gmtime(&time2);
strftime(tbuffer,80,"%Y %B %d %A %I:%M:%S %p UTC",gmt);
seconds = difftime(time2,time1);
fprintf(fptr,"%%========================================================================\n");
fprintf(fptr,"%% end plot \n");
fprintf(fptr,"%% %s\n",tbuffer);
fprintf(fptr,"%% (%.0lf seconds ellapsed)\n",seconds);
fprintf(fptr,"%%========================================================================\n");
fclose(fptr);
return 0;
}
/*-------------------------------------------------------------------------
* plot histogram
*-------------------------------------------------------------------------*/
int plot_ocs_histo(const symseq *x, const long start, const long finish, const time_t time1, const char *plottype, const char *filename, const char *comment,FILE *lptr){
time_t time2;
struct tm *gmt;
double seconds;
int i,mode,Nbin,Nsym;
long n;
char symbol;
char cmdstr[8]="lolli";//"lolli"=default command string
char tbuffer[80];// time buffer
char buf[256];
char *str;
double pos;
const long N=x->getN();
/*----------------------------------------
* parameter data
*----------------------------------------*/
const char plottypelist[][16]={"lolli","lollipop", "die","dice", "spin","spinner", "dna","genome", "dnan","genomen"};
const int Ncmds = sizeof(plottypelist)/sizeof(char)/16;
const int modelist[Ncmds] ={ 0, 0, 1, 1, 2, 2, 3, 3, 4, 4};
const char cmdstrlist[][16] ={"lolli", "die", "spin", "dna", "dnan"};
const int Nmodes = sizeof(cmdstrlist)/sizeof(char)/16;
const char symlist[Nmodes][16] ={"0123456", ".ABCDEF", ".ABCDEF", ".ATCG", ".ATCGN"};
const char yticks[Nmodes][7][16]={
{"0", "1", "2", "3", "4", "5", "6"},
{"0", "\\sfrac{1}{6}", "\\sfrac{1}{3}", "\\sfrac{1}{2}", "\\sfrac{2}{3}", "\\sfrac{5}{6}", "1"},
{"0", "\\sfrac{1}{6}", "\\sfrac{1}{3}", "\\sfrac{1}{2}", "\\sfrac{2}{3}", "\\sfrac{5}{6}", "1"},
{"0", "\\sfrac{1}{4}", "\\sfrac{1}{2}", "\\sfrac{3}{4}", "1", "", "" },
{"0", "\\sfrac{1}{5}", "\\sfrac{2}{5}", "\\sfrac{3}{5}", "\\sfrac{4}{5}", "1", "" },
};
/*----------------------------------------
* process plottype parameter
*----------------------------------------*/
for(i=0,mode=-1;i<Ncmds;i++)
if(strcmp(plottype,plottypelist[i])==0){
mode=modelist[i];
sprintf(cmdstr,cmdstrlist[mode]);
}
if(mode==-1){
mode=modelist[0];
sprintf(cmdstr,cmdstrlist[mode]);
sprintf(buf,"\nERROR: unknown plottype string \"%s\".\n",plottype);
printeft(lptr,buf);
exit(EXIT_FAILURE);
}
Nsym = strlen(symlist[mode]);
for(i=0,Nbin=0,str=(char *)symlist[mode];i<Nsym;i++)if(str[i]!='.')Nbin++;
/*----------------------------------------
* calculate histogram data
*----------------------------------------*/
sprintf(buf,"\nwrite to file %s using plottype \"%s\"; ...\n",filename,cmdstr);printof(lptr,buf);
seqR1 h(Nbin+2);
switch(mode){
case 1: h = ((dieseq *)x)->histogram(start,finish); break;
case 2: h = ((spinseq *)x)->histogram(start,finish); break;
case 3: h = ((dnaseq *)x)->histogram(start,finish); break;
case 4: h = ((dnanseq *)x)->histogram(start,finish); break;
default:
sprintf(buf,"ERROR using function plot_ocs_histo(...): Histogram for the plottype \"%s\" is not supported.\n",plottype);
printeft(lptr,buf);
exit(EXIT_FAILURE);
}
/*----------------------------------------
* generate TeX plot file
*----------------------------------------*/
FILE *fptr;
if(filename==NULL) fptr=stdout;
else fptr=fopen(filename,"w");
if(fptr==NULL){
fprintf(stderr,"plot_die_seq operation failed: Cannot open output file %s .\n",filename);
fprintf(lptr,"plot_die_seq operation failed: Cannot open output file %s .\n",filename);
return -1;
}
if(plot_header(fptr,time1,comment))return -1;
fprintf(fptr,"\\input{shelltop.tex}%%\n");
fprintf(fptr,"\\begin{document}%%\n");
fprintf(fptr," \\gsize%%\n");
fprintf(fptr," \\psset{xunit=10mm,yunit=17mm}%%\n");
fprintf(fptr," \\begin{pspicture}(-0.5,-0.25)(%d,1.16)%%\n",Nbin+1);
fprintf(fptr," \\psset{%%\n");
fprintf(fptr," linecolor=blue,%%\n");
fprintf(fptr," linewidth=0.5pt,%%\n");
fprintf(fptr," plotstyle=bar,%%\n");
fprintf(fptr," barwidth=0.75,%%\n");
fprintf(fptr," fillcolor=cyan,%%\n");
fprintf(fptr," fillstyle=solid,%%\n");
fprintf(fptr," }%%\n");
fprintf(fptr," \\psaxes[linecolor=axis,yAxis=false,showorigin=false,Dx=1,labels=none,ticks=none](0,0)(0,0)(%d,1.16)%% x-axis\n",Nbin+1);
fprintf(fptr," \\psaxes[linecolor=axis,xAxis=false,showorigin=false,Dy=%.4lf,labels=none](0,0)(0,0)(%d,1.16)%% y-axis\n",1.0/(double)Nbin,Nbin+1);
fprintf(fptr," \\psline[linecolor=red,linestyle=dotted,linewidth=1pt](0,%lf)(%d,%lf)%% uniform distribution marker\n",1.0/(double)Nbin,Nbin+1,1.0/(double)Nbin);
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," %% x-axis labelling\n");
fprintf(fptr," %%----------------------------------------\n");
for(i=0;i<Nsym;i++){
symbol = symlist[mode][i];
if(symbol!='.')fprintf(fptr," \\uput{2pt}[-90]{0}(%d,0){$\\%s%c$}%%\n",i,cmdstr,symbol);
}
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," %% histogram data\n");
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," \\savedata{\\pdata}[");
for(n=1;n<=Nbin;n++)fprintf(fptr,"{%ld,%.4lf},",n,h.get(n)/(double)N);
fprintf(fptr,"]%%\n");
fprintf(fptr," \\dataplot{\\pdata}%%\n");
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," %% y-axis labelling\n");
fprintf(fptr," %%----------------------------------------\n");
for(i=0;i<=Nbin;i++){
str=(char *)yticks[mode][i];
//pos=ytickp[mode][i];
pos=(double)i/(double)Nbin;
if(strlen(str)!=0)fprintf(fptr," \\uput[180]{0}(0,%.4lf){$%s$}%%\n",pos,str);
}
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," %% data labelling\n");
fprintf(fptr," %%----------------------------------------\n");
for(n=1;n<=Nbin;n++)fprintf(fptr," \\rput[t](%ld,1.15){$%.0lf$}%%\n",n,h.get(n));
fprintf(fptr," %%\n");
for(n=1;n<=Nbin;n++)fprintf(fptr," \\rput[t](%ld,0.92){$\\scy(%.2lf\\%%)$}%%\n",n,h.get(n)/(double)N*100.0);
fprintf(fptr," \\end{pspicture}%%\n");
fprintf(fptr,"\\end{document}%%\n");
time(&time2); gmt = gmtime(&time2);
strftime(tbuffer,80,"%Y %B %d %A %I:%M:%S %p UTC",gmt);
seconds = difftime(time2,time1);
fprintf(fptr,"%%========================================================================\n");
fprintf(fptr,"%% end data \n");
fprintf(fptr,"%% %s\n",tbuffer);
fprintf(fptr,"%% (%.0lf seconds ellapsed)\n",seconds);
fprintf(fptr,"%%========================================================================\n");
fclose(fptr);
return 0;
}
/*-------------------------------------------------------------------------
* plot auto-correlation for sequence
* Nhw = center lobe half width (used for labelling only)
*-------------------------------------------------------------------------*/
int plot_ocs_auto(const seqR1 *Rxx, const long Nhw, const time_t time1, const char *filename, const char *comment,FILE *lptr){
FILE *fptr;
const long N = Rxx->getN()/2;
const long nSkip=50;
const long Nhres=N/32;//high resolution region = [-Nhres,+Nhres]
const double lobeheight=Rxx->max();
const double psxunit=10,psyunit=10;
const double xscale=3/(double)N;
const double yscale=2.0/lobeheight;
const double u = yscale*Rxx->last();
const double plotmin = (u < -0.15)? u : -0.15;
long n;
int m;
time_t time2;
struct tm *gmt;
double seconds;
char tbuffer[80];
if(filename==NULL) fptr=stdout;
else fptr=fopen(filename,"w");
if(plot_header(fptr,time1,comment))return -1;
fprintf(fptr,"\\input{shelltop.tex}%%\n");
fprintf(fptr,"\\begin{document}%%\n");
fprintf(fptr," \\gsize%%\n");
//fprintf(fptr," \\psset{xunit=3mm,yunit=3mm}%%\n");
fprintf(fptr," \\psset{xunit=%.3lfmm,yunit=%.3lfmm}%%\n",psxunit,psyunit);
fprintf(fptr," \\begin{pspicture}(-3.5,%.3lf)(3.5,2.05)%%\n",plotmin);
fprintf(fptr," %%-------------------------------------\n");
fprintf(fptr," %% parameters\n");
fprintf(fptr," %%-------------------------------------\n");
fprintf(fptr," \\psset{%%\n");
fprintf(fptr," linecolor=blue,%%\n");
fprintf(fptr," linewidth=0.75pt,%%\n");
fprintf(fptr," fillstyle=none,%%\n");
fprintf(fptr," }%%\n");
fprintf(fptr," %%-------------------------------------\n");
fprintf(fptr," %% axes and label positions\n");
fprintf(fptr," %%-------------------------------------\n");
fprintf(fptr," \\psaxes[linecolor=axis,yAxis=false,labels=none,Dx=3]{<->}(0,0)(-3.5,0)(3.5,2.05)%%\n");
fprintf(fptr," \\rput[tl](0.4,2){\\rnode[l]{centerTl}{$\\opair{0}{%.0lf}$}}%% label for top of center lobe\n",Rxx->get(N));
fprintf(fptr," \\rput[r](3.2,0.84){\\rnode[l]{centerRl}{$\\opair{%ld}{%.1lf}$}}%% label for bottom right of center lobe\n",Nhw,Rxx->get(N+Nhw));
fprintf(fptr," \\rput[l](-3.2,0.84){\\rnode[r]{centerLl}{$\\opair{-%ld}{%.1lf}$}}%% label for bottom left of center lobe\n",Nhw,Rxx->get(N-Nhw));
fprintf(fptr," \\rput[l](0.25,0.24){\\rnode[rb]{noiseRl}{$\\opair{%ld}{0}$}}%% label for bottom right of right noise area\n",N);
fprintf(fptr," \\rput[r](-0.25,0.24){\\rnode[lb]{noiseLl}{$\\opair{-%ld}{0}$}}%% label for bottom left of left noise area\n",N);
fprintf(fptr," \\rput[t](-1.7,2){\\rnode[r]{highresl}{\\scs full resolution region}}%% label for full resolution area\n");
fprintf(fptr," \\rput[b](-1.6,1.2){\\rnode[b]{lowresLl}{\\scs1:%ld resolution region}}%% label for left low resolution area\n",nSkip);
fprintf(fptr," \\rput[b](1.6,1.2){\\rnode[b]{lowresRl}{\\scs1:%ld resolution region}}%% label for right low resolution area\n",nSkip);
fprintf(fptr," %%-------------------------------------\n");
fprintf(fptr," %% node definitions\n");
fprintf(fptr," %%-------------------------------------\n");
fprintf(fptr," \\psset{xunit=%.10lfmm,yunit=%.10lfmm}%%\n",xscale*psxunit,yscale*psyunit);
fprintf(fptr," \\pnode(0,%.4lf){centerTp}%% top of center lobe\n",lobeheight);
fprintf(fptr," \\pnode(%ld,%.4lf){centerRp}%% bottom right of center lobe\n",Nhw,Rxx->get(N+Nhw));
fprintf(fptr," \\pnode(-%ld,%.4lf){centerLp}%% bottom left of center lobe\n",Nhw,Rxx->get(N-Nhw));
fprintf(fptr," \\pnode(%ld,0){noiseRp}%% bottom right of noise lobe\n",N);
fprintf(fptr," \\pnode(-%ld,0){noiseLp}%% bottom left of noise lobe\n",N);
fprintf(fptr," \\pnode(-%ld,%.4lf){highresp}%% high resolution left boundary\n",Nhres,lobeheight*0.85);
fprintf(fptr," %%-------------------------------------\n");
fprintf(fptr," %% guide markings\n");
fprintf(fptr," %%-------------------------------------\n");
fprintf(fptr," \\psline[linestyle=dotted,linecolor=red](-%ld,0)(-%ld,%.3lf)%%\n",N,N,lobeheight);
fprintf(fptr," \\psline[linestyle=dotted,linecolor=red](+%ld,0)(+%ld,%.3lf)%%\n",N,N,lobeheight);
fprintf(fptr," \\psline[linestyle=dotted,linecolor=red](-%ld,0)(-%ld,%.3lf)%%\n",Nhres,Nhres,lobeheight);
fprintf(fptr," \\psline[linestyle=dotted,linecolor=red](+%ld,0)(+%ld,%.3lf)%%\n",Nhres,Nhres,lobeheight);
fprintf(fptr," \\ncline[linecolor=red]{->}{centerTl}{centerTp}%%\n");
fprintf(fptr," \\ncline[linecolor=red]{->}{centerRl}{centerRp}%%\n");
fprintf(fptr," \\ncline[linecolor=red]{->}{centerLl}{centerLp}%%\n");
fprintf(fptr," \\ncarc[linecolor=red,arcangle=-22]{->}{noiseRl}{noiseRp}%%\n");
fprintf(fptr," \\ncarc[linecolor=red,arcangle=22]{->}{noiseLl}{noiseLp}%%\n");
fprintf(fptr," \\ncline[linecolor=red]{->}{highresl}{highresp}%%\n");
fprintf(fptr," %%-----------------------------------------\n");
fprintf(fptr," %% plot left side data at 1:%ld resolution\n",nSkip);
fprintf(fptr," %%-----------------------------------------\n");
fprintf(fptr," \\psline(-%.3lf,%.3lf)",3.4/xscale,Rxx->first());//left side outlying region
for(n=0,m=0;n<=(N-Nhres);n+=nSkip,m++){
if(m%5==0)fprintf(fptr,"%%\n ");
fprintf(fptr,"(%5ld,%10.3lf)",n-N,Rxx->get(n));
}
fprintf(fptr," (%5ld,%10.3lf)%%\n",-Nhres,Rxx->get(N-Nhres));
fprintf(fptr," %%-----------------------------------------\n");
fprintf(fptr," %% plot center data at full resolution\n");
fprintf(fptr," %%-----------------------------------------\n");
fprintf(fptr," \\psline");
for(n=(N-Nhres),m=0;n<=(N+Nhres);n++,m++){
if(m%5==0)fprintf(fptr,"%%\n ");
fprintf(fptr,"(%5ld,%10.3lf)",n-N,Rxx->get(n));
}
fprintf(fptr,"%%\n");
fprintf(fptr," %%-----------------------------------------\n");
fprintf(fptr," %% plot right side data at 1:%ld resolution\n",nSkip);
fprintf(fptr," %% (plotted in reverse for symmetry with left side)\n");
fprintf(fptr," %%-----------------------------------------\n");
fprintf(fptr," \\psline(%.3lf,%.3lf)",3.4/xscale,Rxx->last());//right side outlying region
for(n=Rxx->getN()-1,m=0;n>=N+Nhres;n-=nSkip,m++){
if(m%5==0)fprintf(fptr,"%%\n ");
fprintf(fptr,"(%5ld,%10.3lf)",n-N,Rxx->get(n));
}
fprintf(fptr,"(%5ld,%10.3lf)%%\n",Nhres,Rxx->get(N+Nhres));
fprintf(fptr," %%-----------------------------------------\n");
fprintf(fptr," %% re-plot left side of center spike for better clarity\n");
fprintf(fptr," %%-----------------------------------------\n");
fprintf(fptr," \\psline(-1,%10.3lf)(0,%10.3lf)\n",Rxx->get(N-1),Rxx->get(N));
fprintf(fptr," \\end{pspicture}%%\n");
fprintf(fptr,"\\end{document}%%\n");
time(&time2);
gmt = gmtime(&time2);
strftime(tbuffer,80,"%Y %B %d %A %I:%M:%S %p UTC",gmt);
seconds = difftime(time2,time1);
fprintf(fptr,"%%========================================================================\n");
fprintf(fptr,"%% end data \n");
fprintf(fptr,"%% %s\n",tbuffer);
fprintf(fptr,"%% (%.0lf seconds ellapsed)\n",seconds);
fprintf(fptr,"%%========================================================================\n");
fclose(fptr);
return 0;
}
/*-------------------------------------------------------------------------
* plot a sequence <x> from <start> to <finish>
* and downsampled by a factor of <M>
* with the y-axis ranging from <min> to <max>
* and plot size measuring <widthmm> mm by <heightmm> mm
* and write the LaTeX plot file to <filename>
* with <comment> incluced as a comment.
*-------------------------------------------------------------------------*/
int plot_R1_seq(const seqR1 *x, const long start, const long finish, const long M, const double min, const double max, const double widthmm, const double heightmm, const char *filename, const char *comment){
time_t time1;time(&time1); //starting time stamp
const long N=x->getN();
const double xunit=(widthmm-20)/((double)finish); //(mm/unit)
const double yunit=(heightmm-4)/(max-min); //(mm/unit)
const char sublinestyle[16]="dotted";
const long xgridstep = (long)pow(10.0,floor(log10((double)(finish+1))+0.5)-1);
const long xsubgridstep = xgridstep/10;
long xlabelstep = xgridstep;
const double ygridstep = pow(10.0,floor(log10(max-min))-1);
const double axxb=(double)(finish+1); //axes rectangle from (axxa,axya) to (axxb,axyb)
const double axxa=-10.0/xunit;
const double axyb=max+ygridstep;
const double axya=min-ygridstep;
const double picxb=(double)finish+8.0/xunit; //right side of picture is 10mm past last data point
const double picxa=-12.0/xunit;
const double picyb=axyb*1.05;
const double picya=axya*1.4;
long n,m;
double xn,y;
FILE *fptr;
// const long M=1;
/*----------------------------------------
* open plot file
*----------------------------------------*/
fptr=plot_open(filename,time1,comment);
/*----------------------------------------
* generate TeX plot file
*----------------------------------------*/
fprintf(fptr,"\\input{shelltop.tex}%%\n");
fprintf(fptr,"\\begin{document}%%\n");
fprintf(fptr," \\gsize%%\n");
fprintf(fptr," \\psset{labelsep=1pt,linewidth=0.75pt,xunit=%.6lfmm,yunit=%.6lfmm}%%\n",xunit,yunit);
fprintf(fptr," \\begin{pspicture}(%.3lf,%.3lf)(%.3lf,%.3lf)%%\n",picxa,picya,picxb,picyb);
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," %% axes\n");
fprintf(fptr," %%----------------------------------------\n");
//fprintf(fptr," \\psaxes[linecolor=axis,yAxis=false,showorigin=false,Dx=100]{->}(0,0)(%.3lf,%.3lf)(%.3lf,%.3lf)%% x-axis\n",axxa,axya,axxb,axyb);
if(min<0)fprintf(fptr," \\psaxes[linecolor=axis,xAxis=false,showorigin=false,labels=none]{<->}(0,0)(%.3lf,%.3lf)(%.3lf,%.3lf)%% y-axis\n",axxa,axya,axxb,axyb);
else fprintf(fptr," \\psaxes[linecolor=axis,xAxis=false,showorigin=false,labels=none]{->}(0,0)(%.3lf,0)(%.3lf,%.3lf)%% y-axis\n",axxa,axxb,axyb);
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," %% horizontal graph lines\n");
fprintf(fptr," %%----------------------------------------\n");
for(y=max,n=0;y>=min;y-=(max/10),n++)
if(n%10==0)fprintf(fptr," \\psline[linecolor=graph,linewidth=1.0pt,linestyle=solid](0,%.3lf)(%.3lf,%.3lf)%%\n",y,axxb,y);
else fprintf(fptr," \\psline[linecolor=graph,linewidth=0.5pt,linestyle=%s](0,%.3lf)(%.3lf,%.3lf)%%\n",sublinestyle,y,axxb,y);
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," %% y-axis labels\n");
fprintf(fptr," %%----------------------------------------\n");
for(y=max;y>=min;y-=(max/5))
if(y==0)fprintf(fptr," \\uput{3pt}[180]{0}(0,%.3lf){0}%%\n",y);
else fprintf(fptr," \\uput{3pt}[180]{0}(0,%.3lf){%+.3lf}%%\n",y,y);
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," %% vertical graph lines\n");
fprintf(fptr," %%----------------------------------------\n");
for(n=xsubgridstep;n<=(finish+1);n+=xsubgridstep)
if(n%xgridstep==0)fprintf(fptr," \\psline[linecolor=graph,linewidth=1.0pt,linestyle=%-6s](%6ld,%.3lf)(%6ld,%.3lf)%%\n","solid",n,min,n,max);
else fprintf(fptr," \\psline[linecolor=graph,linewidth=0.5pt,linestyle=%-6s](%6ld,%.3lf)(%6ld,%.3lf)%%\n",sublinestyle,n,min,n,max);
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," %% labels along top of graph\n");
fprintf(fptr," %%----------------------------------------\n");
if(finish/xlabelstep>10)xlabelstep*=2;
if(finish/xlabelstep>10)xlabelstep*=2;
if(finish/xlabelstep>10)xlabelstep*=2;
for(n=xlabelstep;n<=(finish+1);n+=xlabelstep)
fprintf(fptr," \\rput[t](%ld,%.3lf){%ld}%%\n",n,axyb,n);
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," %% data points\n");
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," \\psline");
for(n=start,m=0;n<=finish;n+=M,m++){
xn=x->get(n);
fprintf(fptr,"(%5ld,%12.6lf)",n,xn);
if(m%5==4&&n!=finish)fprintf(fptr,"%%\n ");
}
fprintf(fptr,"\n");
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," %% fractional labels along bottom of plot\n");
fprintf(fptr," %%----------------------------------------\n");
for(n=N/10,m=1;n<=(finish+1);n+=N/10,m++){
if (n==N/10)fprintf(fptr," \\rput[b](%5ld,%.3lf){$\\frac{ \\xN}{10}$}%%\n",n,axya);
else if(m==10) fprintf(fptr," \\rput[b](%5ld,%.3lf){$\\xN$}%%\n",N,axya);
else fprintf(fptr," \\rput[b](%5ld,%.3lf){$\\frac{%2ld\\xN}{10}$}%%\n",n,axya,m);
}
fprintf(fptr," \\end{pspicture}%%\n");
fprintf(fptr,"\\end{document}%%\n");
/*----------------------------------------
* close plot file
*----------------------------------------*/
plot_close(fptr,time1);
return 0;
}
/*-------------------------------------------------------------------------
* plot DFT sequence
* with full resolutoin in the +/- 50 sample area around hres1, hres2, hres3,
* and with other areas downsampled by a factor of M.
* If M<2, then use full resolution everywhere.
*-------------------------------------------------------------------------*/
int plot_dft_seq(const seqR1 *x, const long start, const long finish,const double mindB, const double maxdB, const long M, const long hres1, const long hres2, const long hres3, const char *filename, const char *comment){
time_t time1;time(&time1); //starting time stamp
const long N=x->getN();
const double xunit=140.0/((double)finish); //(mm/unit)
const double picxb=(double)finish+8.0/xunit; //right side of picture is 10mm past last data point
const double picxa=-10.0/xunit;
const double picyb=maxdB+1+0.1;
const double picya=mindB-1-0.1;
const double axxb=(double)finish; //axes rectangle from (axxa,axya) to (axxb,axyb)
const double axxa=-10.0/xunit;
const double axyb=maxdB+1;
const double axya=mindB-1;
long n,m;
double xn;
FILE *fptr;
double y;
const char sublinestyle[16]="dotted";
const long xgridstep = (long)pow(10.0,floor(log10((double)(finish+1))+0.5)-1);
const long xsubgridstep = xgridstep/10;
long xlabelstep = xgridstep;
//const double ygridstep = pow(10.0,floor(log10(maxdB-mindB))-1);
/*----------------------------------------
* open plot file
*----------------------------------------*/
fptr=plot_open(filename,time1,comment);
/*----------------------------------------
* generate TeX plot file
*----------------------------------------*/
fprintf(fptr,"\\input{shelltop.tex}%%\n");
fprintf(fptr,"\\begin{document}%%\n");
fprintf(fptr," \\gsize%%\n");
fprintf(fptr," \\psset{labelsep=1pt,linewidth=0.75pt,xunit=%.3lfmm,yunit=3mm}%%\n",xunit);
fprintf(fptr," \\begin{pspicture}(%.3lf,%.3lf)(%.3lf,%.3lf)%%\n",picxa,picya,picxb,picyb);
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," %% axes\n");
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," \\psaxes[linecolor=axis,xAxis=false,showorigin=false,labels=none]{<->}(0,0)(%.3lf,%.3lf)(%.3lf,%.3lf)%% y-axis\n",axxa,axya,axxb,axyb);
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," %% horizontal graph lines\n");
fprintf(fptr," %%----------------------------------------\n");
for(n=(long)maxdB;n>=(long)mindB;n--)fprintf(fptr," \\psline[linecolor=graph,linewidth=0.5pt,linestyle=dotted](0,%2ld)(%ld,%2ld)%%\n",n,finish,n);
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," %% y-axis dB labels\n");
fprintf(fptr," %%----------------------------------------\n");
for(n=(long)maxdB;n>=(long)mindB;n--)fprintf(fptr," \\uput{3pt}[180]{0}(0,%2ld){%+2lddB}%%\n",n,n);
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," %% vertical graph lines\n");
fprintf(fptr," %%----------------------------------------\n");
for(n=xsubgridstep;n<=(finish+1);n+=xsubgridstep)
if(n%xgridstep==0)fprintf(fptr," \\psline[linecolor=graph,linewidth=1.0pt,linestyle=%-6s](%6ld,%.3lf)(%6ld,%.3lf)%%\n","solid",n,mindB,n,maxdB);
else fprintf(fptr," \\psline[linecolor=graph,linewidth=0.5pt,linestyle=%-6s](%6ld,%.3lf)(%6ld,%.3lf)%%\n",sublinestyle,n,mindB,n,maxdB);
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," %% Hz labels along top of graph\n");
fprintf(fptr," %%----------------------------------------\n");
if(finish/xlabelstep>10)xlabelstep*=2;
if(finish/xlabelstep>10)xlabelstep*=2;
if(finish/xlabelstep>10)xlabelstep*=2;
for(n=xlabelstep;n<=(finish+1);n+=xlabelstep)
fprintf(fptr," \\rput[t](%ld,%.3lf){%ld{\\scriptsize Hz}}%%\n",n,axyb,n);
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," %% data points\n");
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," \\psline");
for(n=start,m=0;n<=finish;m++){
xn=x->get(n);
fprintf(fptr,"(%5ld,%12.6lf)",n,10.0*log10(xn));
if(m%5==4&&n!=finish)fprintf(fptr,"%%\n ");
if (M<2) n++; // If M<2, use full resolution everywhere
else if((n>=hres1-50)&&(n<=hres1+50)) n++; // use full resolution in the area around hres1
else if((n>=hres2-50)&&(n<=hres2+50)) n++; // use full resolution in the area around hres2
else if((n>=hres3-50)&&(n<=hres3+50)) n++; // use full resolution in the area around hres3
else n+=M; // else downsample by factor M
}
fprintf(fptr,"%%\n");
fprintf(fptr," \\uput[0]{0}(%5ld,0){{\\Large\\color{blue}$\\cdots$}}%%\n",n);
fprintf(fptr," %%----------------------------------------\n");
fprintf(fptr," %% pi labels along bottom of plot\n");
fprintf(fptr," %%----------------------------------------\n");
for(y=(double)N/20.0,n=1;y<=(double)(finish+2);y+=(double)N/20.0,n++){
if (n==1) fprintf(fptr," \\rput[b](%10.3lf,%+8.3lf){$\\frac{ \\pi}{10}$}%%\n",y,axya);
else if(n==10)fprintf(fptr," \\rput[b](%10.3lf,%+8.3lf){$\\pi$}%%\n",y,axya);
else if(n==20)fprintf(fptr," \\rput[b](%10.3lf,%+8.3lf){$2\\pi$}%%\n",y,axya);
else fprintf(fptr," \\rput[b](%10.3lf,%+8.3lf){$\\frac{%ld\\pi}{10}$}%%\n",y,axya,n);
}
fprintf(fptr," \\end{pspicture}%%\n");
fprintf(fptr,"\\end{document}%%\n");
/*----------------------------------------
* close plot file
*----------------------------------------*/
plot_close(fptr,time1);
return 0;
}
/*-------------------------------------------------------------------------
* printf string <*str> to standard output and to file <*fptr>
* write a time stamp into file as well
*-------------------------------------------------------------------------*/
void printeft(FILE *fptr, const char *str){
time_t time2;time(&time2);
struct tm *gmt;
gmt = gmtime(&time2);
char tbuffer[80];
strftime(tbuffer,80,"%Y %B %d %A %I:%M:%S %p UTC",gmt);
fprintf(stderr,"%s [%s]",str,tbuffer);
fprintf(fptr,"%s [%s]",str,tbuffer);
}
/*-------------------------------------------------------------------------
* printf string <*str> to standard output and to file <*fptr>
* write a time stamp into file as well
*-------------------------------------------------------------------------*/
void printoft(FILE *fptr, const char *str){
time_t time2;time(&time2);
struct tm *gmt;
gmt = gmtime(&time2);
char tbuffer[80];
strftime(tbuffer,80,"%Y %B %d %A %I:%M:%S %p UTC",gmt);
fprintf(stdout,"%s [%s]",str,tbuffer);
fprintf(fptr,"%s [%s]",str,tbuffer);
}
/*-------------------------------------------------------------------------
* printf string <*str> time stamp
* to standard output and to file <*fptr>
*-------------------------------------------------------------------------*/
void printoft(FILE *fptr, const char *str, const time_t time1){
time_t time2;time(&time2);
struct tm *gmt;
gmt = gmtime(&time2);
char tbuffer[80];
strftime(tbuffer,80,"%Y %B %d %A %I:%M:%S %p UTC",gmt);
fprintf(stdout,"%s [%s]\n",str,tbuffer);
fprintf(fptr,"%s [%s]\n",str,tbuffer);
}
/*-------------------------------------------------------------------------
* printf string <*str> and time ellapsed record
* to standard output and to file <*fptr>
*-------------------------------------------------------------------------*/
void printofe(FILE *fptr, const char *str, const time_t time1){
time_t time2;time(&time2);
double seconds=difftime(time2,time1);
struct tm *gmt;
gmt = gmtime(&time2);
char tbuffer[80];
strftime(tbuffer,80,"%Y %B %d %A %I:%M:%S %p UTC",gmt);
fprintf(stdout,"%s [%.0lf seconds elapsed]\n",str,seconds);
fprintf(fptr,"%s [%.0lf seconds elapsed]\n",str,seconds);
}