comedi_calibrate/comedi_calibrate.c

   1 /*
   2    A little auto-calibration utility, for boards
   3    that support it.
   4
   5    Right now, it only supports NI E series boards,
   6    but it should be easily portable.
   7
   8    A few things need improvement here:
   9     - current system gets "close", but doesn't
  10       do any fine-tuning
  11     - no pre/post gain discrimination for the
  12       A/D zero offset.
  13     - should read (and use) the actual reference
  14       voltage value from eeprom
  15     - statistics would be nice, to show how good
  16       the calibration is.
  17     - doesn't check unipolar ranges
  18     - "alternate calibrations" would be cool--to
  19       accurately measure 0 in a unipolar range
  20     - more portable
  21  */
  22
  23 #define _GNU_SOURCE
  24
  25 #include <stdio.h>
  26 #include <comedilib.h>
  27 #include <fcntl.h>
  28 #include <unistd.h>
  29 #include <errno.h>
  30 #include <getopt.h>
  31 #include <ctype.h>
  32 #include <math.h>
  33 #include <stdlib.h>
  34 #include <string.h>
  35
  36 #include "calib.h"
  37
  38 /* global variables */
  39
  40 caldac caldacs[N_CALDACS];
  41 int n_caldacs;
  42
  43 observable observables[N_OBSERVABLES];
  44 int n_observables;
  45
  46 comedi_t *dev;
  47
  48 int ad_subdev;
  49 int da_subdev;
  50 int eeprom_subdev;
  51 int caldac_subdev;
  52
  53 char *drivername = NULL;
  54 char *devicename = NULL;
  55
  56 int verbose = 0;
  57
  58 int device_status = STATUS_UNKNOWN;
  59
  60 /* */
  61
  62
  63 struct board_struct{
  64         char *name;
  65         void (*setup)(void);
  66 };
  67
  68 struct board_struct drivers[] = {
  69         { "ni_pcimio",          ni_setup },
  70         { "ni_atmio",           ni_setup },
  71         { "ni_mio_cs",          ni_setup },
  72 };
  73 #define n_drivers (sizeof(drivers)/sizeof(drivers[0]))
  74
  75 int do_dump = 0;
  76 int do_reset = 1;
  77 int do_calibrate = 1;
  78 int do_results = 0;
  79 int do_output = 1;
  80
  81 struct option options[] = {
  82         { "verbose", 0, 0, 'v' },
  83         { "quiet", 0, 0, 'q' },
  84         { "file", 1, 0, 'f' },
  85         { "driver-name", 1, 0, 0x1000 },
  86         { "device-name", 1, 0, 0x1001 },
  87         { "reset", 0, &do_reset, 1 },
  88         { "no-reset", 0, &do_reset, 0 },
  89         { "calibrate", 0, &do_calibrate, 1 },
  90         { "no-calibrate", 0, &do_calibrate, 0 },
  91         { "dump", 0, &do_dump, 1 },
  92         { "no-dump", 0, &do_dump, 0 },
  93         { "results", 0, &do_results, 1 },
  94         { "no-results", 0, &do_results, 0 },
  95         { "output", 0, &do_output, 1 },
  96         { "no-output", 0, &do_output, 0 },
  97         { 0 },
  98 };
  99
 100 int main(int argc, char *argv[])
 101 {
 102         char *fn = NULL;
 103         int c;
 104         int i;
 105         struct board_struct *this_board;
 106         int index;
 107
 108         fn = "/dev/comedi0";
 109         while (1) {
 110                 c = getopt_long(argc, argv, "f:vq", options, &index);
 111                 if (c == -1)break;
 112                 switch (c) {
 113                 case 0:
 114                         continue;
 115                 case 'f':
 116                         fn = optarg;
 117                         break;
 118                 case 'v':
 119                         verbose++;
 120                         break;
 121                 case 'q':
 122                         verbose--;
 123                         break;
 124                 case 0x1000:
 125                         drivername = optarg;
 126                         break;
 127                 case 0x1001:
 128                         devicename = optarg;
 129                         break;
 130                 default:
 131                         printf("bad option %d\n",c);
 132                         exit(1);
 133                 }
 134         }
 135
 136         dev = comedi_open(fn);
 137         if (dev == NULL ) {
 138                 perror(fn);
 139                 exit(0);
 140         }
 141
 142         if(!drivername)
 143                 drivername=comedi_get_driver_name(dev);
 144         if(!devicename)
 145                 devicename=comedi_get_board_name(dev);
 146
 147         ad_subdev=comedi_find_subdevice_by_type(dev,COMEDI_SUBD_AI,0);
 148         da_subdev=comedi_find_subdevice_by_type(dev,COMEDI_SUBD_AO,0);
 149         caldac_subdev=comedi_find_subdevice_by_type(dev,COMEDI_SUBD_CALIB,0);
 150         eeprom_subdev=comedi_find_subdevice_by_type(dev,COMEDI_SUBD_MEMORY,0);
 151
 152         for(i=0;i<n_drivers;i++){
 153                 if(!strcmp(drivers[i].name,drivername)){
 154                         this_board = drivers+i;
 155                         goto ok;
 156                 }
 157         }
 158         printf("Driver %s unknown\n",drivername);
 159         return 1;
 160
 161 ok:
 162         this_board->setup();
 163
 164         if(device_status<STATUS_DONE){
 165                 printf("Warning: device not fully calibrated due to insufficient information\n");
 166                 printf("Please send this output to <ds@schleef.org>\n");
 167                 if(verbose<1)verbose=1;
 168                 if(device_status==STATUS_UNKNOWN){
 169                         do_reset=1;
 170                         do_dump=1;
 171                         do_calibrate=0;
 172                         do_results=0;
 173                 }
 174                 if(device_status==STATUS_SOME){
 175                         do_reset=1;
 176                         do_dump=1;
 177                         do_calibrate=1;
 178                         do_results=1;
 179                 }
 180                 if(device_status==STATUS_GUESS){
 181                         do_reset=1;
 182                         do_dump=1;
 183                         do_calibrate=1;
 184                         do_results=1;
 185                 }
 186         }
 187         if(verbose>=0){
 188                 printf("$Id$\n");
 189                 printf("Driver name: %s\n",drivername);
 190                 printf("Device name: %s\n",devicename);
 191                 printf("Comedi version: %d.%d.%d\n",
 192                         (comedi_get_version_code(dev)>>16)&0xff,
 193                         (comedi_get_version_code(dev)>>8)&0xff,
 194                         (comedi_get_version_code(dev))&0xff);
 195         }
 196
 197         if(do_reset)reset_caldacs();
 198         if(do_dump)observe();
 199         if(do_calibrate && do_cal)do_cal();
 200         if(do_results)observe();
 201
 202         return 0;
 203 }
 204
 205 void set_target(int obs,double target)
 206 {
 207         comedi_range *range;
 208         lsampl_t maxdata, data;
 209
 210         range = comedi_get_range(dev,observables[obs].preobserve_insn.subdev,
 211                 CR_CHAN(observables[obs].preobserve_insn.chanspec),
 212                 CR_RANGE(observables[obs].preobserve_insn.chanspec));
 213         maxdata = comedi_get_maxdata(dev,observables[obs].preobserve_insn.subdev,
 214                 CR_CHAN(observables[obs].preobserve_insn.chanspec));
 215
 216         data = comedi_from_phys(target,range,maxdata);
 217
 218         observables[obs].preobserve_data = data;
 219         observables[obs].target = comedi_to_phys(data,range,maxdata);
 220 }
 221
 222 void observe(void)
 223 {
 224         int i;
 225
 226         for(i=0;i<n_observables;i++){
 227                 preobserve(i);
 228                 DPRINT(0,"%s\n",observables[i].name);
 229                 measure_observable(i);
 230                 if(verbose>=1){
 231                         observable_dependence(i);
 232                 }
 233         }
 234
 235 }
 236
 237 void preobserve(int obs)
 238 {
 239         if(observables[obs].preobserve_insn.n!=0){
 240                 comedi_do_insn(dev,&observables[obs].preobserve_insn);
 241         }
 242 }
 243
 244 void measure_observable(int obs)
 245 {
 246         char s[32];
 247         int n;
 248         new_sv_t sv;
 249
 250         new_sv_init(&sv,dev,
 251                 observables[obs].observe_insn.subdev,
 252                 CR_CHAN(observables[obs].observe_insn.chanspec),
 253                 CR_RANGE(observables[obs].observe_insn.chanspec),
 254                 CR_AREF(observables[obs].observe_insn.chanspec));
 255         sv.order=7;
 256         n=new_sv_measure(&sv);
 257
 258         sci_sprint_alt(s,sv.average,sv.error);
 259         DPRINT(0,"offset %s, target %g\n",s,observables[obs].target);
 260 }
 261
 262 void observable_dependence(int obs)
 263 {
 264         int i;
 265         linear_fit_t l;
 266
 267         for(i=0;i<n_caldacs;i++){
 268                 check_gain_chan_x(&l,
 269                         observables[obs].observe_insn.chanspec,i);
 270         }
 271
 272 }
 273
 274
 275 void postgain_cal(int obs1, int obs2, int dac)
 276 {
 277         double offset1,offset2;
 278         linear_fit_t l;
 279         double slope1,slope2;
 280         double a;
 281         double gain;
 282         comedi_range *range1,*range2;
 283
 284         DPRINT(0,"postgain: %s; %s\n",observables[obs1].name,
 285                 observables[obs2].name);
 286         preobserve(obs1);
 287         check_gain_chan_x(&l,observables[obs1].observe_insn.chanspec,dac);
 288         offset1=linear_fit_func_y(&l,caldacs[dac].current);
 289         DPRINT(2,"obs1: [%d] offset %g\n",obs1,offset1);
 290         range1 = comedi_get_range(dev,observables[obs1].observe_insn.subdev,
 291                 CR_CHAN(observables[obs1].observe_insn.chanspec),
 292                 CR_RANGE(observables[obs1].observe_insn.chanspec));
 293         slope1=l.slope;
 294
 295         preobserve(obs2);
 296         check_gain_chan_x(&l,observables[obs2].observe_insn.chanspec,dac);
 297         offset2=linear_fit_func_y(&l,caldacs[dac].current);
 298         DPRINT(2,"obs2: [%d] offset %g\n",obs2,offset2);
 299         range2 = comedi_get_range(dev,observables[obs2].observe_insn.subdev,
 300                 CR_CHAN(observables[obs2].observe_insn.chanspec),
 301                 CR_RANGE(observables[obs2].observe_insn.chanspec));
 302         slope2=l.slope;
 303
 304         gain = (range1->max-range1->min)/(range2->max-range2->min);
 305         DPRINT(4,"range1 %g range2 %g\n", range1->max-range1->min,
 306                 range2->max-range2->min);
 307         DPRINT(3,"gain: %g\n",gain);
 308
 309         DPRINT(3,"difference: %g\n",offset2-offset1);
 310
 311         a = (offset1-offset2)/(slope1-slope2);
 312         a=caldacs[dac].current-a;
 313
 314         caldacs[dac].current=rint(a);
 315         update_caldac(dac);
 316         usleep(100000);
 317
 318         DPRINT(0,"caldac[%d] set to %g (%g)\n",dac,rint(a),a);
 319
 320         if(verbose>=2){
 321                 preobserve(obs1);
 322                 measure_observable(obs1);
 323                 preobserve(obs2);
 324                 measure_observable(obs2);
 325         }
 326 }
 327
 328 void cal1(int obs, int dac)
 329 {
 330         linear_fit_t l;
 331         double a;
 332
 333         DPRINT(0,"linear: %s\n",observables[obs].name);
 334         preobserve(obs);
 335         check_gain_chan_x(&l,observables[obs].observe_insn.chanspec,dac);
 336         a=linear_fit_func_x(&l,observables[obs].target);
 337
 338         caldacs[dac].current=rint(a);
 339         update_caldac(dac);
 340         usleep(100000);
 341
 342         DPRINT(0,"caldac[%d] set to %g (%g)\n",dac,rint(a),a);
 343         if(verbose>=3){
 344                 measure_observable(obs);
 345         }
 346 }
 347
 348 void cal1_fine(int obs, int dac)
 349 {
 350         linear_fit_t l;
 351         double a;
 352
 353         DPRINT(0,"linear fine: %s\n",observables[obs].name);
 354         preobserve(obs);
 355         check_gain_chan_fine(&l,observables[obs].observe_insn.chanspec,dac);
 356         a=linear_fit_func_x(&l,observables[obs].target);
 357
 358         caldacs[dac].current=rint(a);
 359         update_caldac(dac);
 360         usleep(100000);
 361
 362         DPRINT(0,"caldac[%d] set to %g (%g)\n",dac,rint(a),a);
 363         if(verbose>=3){
 364                 measure_observable(obs);
 365         }
 366 }
 367
 368 #if 0
 369 void chan_cal(int adc,int cdac,int range,double target)
 370 {
 371         linear_fit_t l;
 372         double offset;
 373         double gain;
 374         double a;
 375         char s[32];
 376
 377         check_gain_chan_x(&l,CR_PACK(adc,range,AREF_OTHER),cdac);
 378         offset=linear_fit_func_y(&l,caldacs[cdac].current);
 379         gain=l.slope;
 380
 381         a=caldacs[cdac].current+(target-offset)/gain;
 382
 383         caldacs[cdac].current=rint(a);
 384         update_caldac(cdac);
 385
 386         read_chan2(s,adc,range);
 387         DPRINT(1,"caldac[%d] set to %g, offset=%s\n",cdac,a,s);
 388 }
 389 #endif
 390
 391
 392 #if 0
 393 void channel_dependence(int adc,int range)
 394 {
 395         int i;
 396         double gain;
 397
 398         for(i=0;i<n_caldacs;i++){
 399                 gain=check_gain_chan(adc,range,i);
 400         }
 401 }
 402 #endif
 403
 404 #if 0
 405 void caldac_dependence(int caldac)
 406 {
 407         int i;
 408         double gain;
 409
 410         for(i=0;i<8;i++){
 411                 gain=check_gain_chan(i,0,caldac);
 412         }
 413 }
 414 #endif
 415
 416
 417 void setup_caldacs(void)
 418 {
 419         int n_chan,i;
 420
 421         if(caldac_subdev<0){
 422                 printf("no calibration subdevice\n");
 423                 return;
 424         }
 425
 426         n_chan=comedi_get_n_channels(dev,caldac_subdev);
 427
 428         for(i=0;i<n_chan;i++){
 429                 caldacs[i].subdev=caldac_subdev;
 430                 caldacs[i].chan=i;
 431                 caldacs[i].maxdata=comedi_get_maxdata(dev,caldac_subdev,i);
 432                 caldacs[i].current=0;
 433         }
 434
 435         n_caldacs=n_chan;
 436 }
 437
 438 void reset_caldacs(void)
 439 {
 440         int i;
 441
 442         for(i=0;i<n_caldacs;i++){
 443                 caldacs[i].current=caldacs[i].maxdata/2;
 444                 update_caldac(i);
 445         }
 446 }
 447
 448 void update_caldac(int i)
 449 {
 450         int ret;
 451
 452         DPRINT(4,"update %d %d %d\n",caldacs[i].subdev,caldacs[i].chan,caldacs[i].current);
 453         if(caldacs[i].current<0){
 454                 DPRINT(1,"caldac set out of range (%d<0)\n",caldacs[i].current);
 455                 caldacs[i].current=0;
 456         }
 457         if(caldacs[i].current>caldacs[i].maxdata){
 458                 DPRINT(1,"caldac set out of range (%d>%d)\n",
 459                         caldacs[i].current,caldacs[i].maxdata);
 460                 caldacs[i].current=caldacs[i].maxdata;
 461         }
 462
 463         ret = comedi_data_write(dev,caldacs[i].subdev,caldacs[i].chan,0,0,
 464                 caldacs[i].current);
 465         if(ret<0)perror("update_caldac()");
 466 }
 467
 468 #if 0
 469 void check_gain(int ad_chan,int range)
 470 {
 471         int i;
 472
 473         for(i=0;i<n_caldacs;i++){
 474                 check_gain_chan(ad_chan,range,i);
 475         }
 476 }
 477 #endif
 478
 479 #if 0
 480 double check_gain_chan(int ad_chan,int range,int cdac)
 481 {
 482         linear_fit_t l;
 483
 484         return check_gain_chan_x(&l,CR_PACK(ad_chan,range,AREF_OTHER),cdac);
 485 }
 486 #endif
 487
 488
 489 double check_gain_chan_x(linear_fit_t *l,unsigned int ad_chanspec,int cdac)
 490 {
 491         int orig,i,n;
 492         int step;
 493         new_sv_t sv;
 494         double sum_err;
 495         int sum_err_count=0;
 496         char str[20];
 497
 498         n=caldacs[cdac].maxdata+1;
 499         memset(l,0,sizeof(*l));
 500
 501         step=n/256;
 502         if(step<1)step=1;
 503         l->n=0;
 504
 505         l->y_data=malloc(n*sizeof(double)/step);
 506         if(l->y_data == NULL)
 507         {
 508                 perror("comedi_calibrate");
 509                 exit(1);
 510         }
 511
 512         orig=caldacs[cdac].current;
 513
 514         new_sv_init(&sv,dev,0,
 515                 CR_CHAN(ad_chanspec),
 516                 CR_RANGE(ad_chanspec),
 517                 CR_AREF(ad_chanspec));
 518
 519         caldacs[cdac].current=0;
 520         update_caldac(cdac);
 521         usleep(100000);
 522
 523         new_sv_measure(&sv);
 524
 525         sum_err=0;
 526         for(i=0;i*step<n;i++){
 527                 caldacs[cdac].current=i*step;
 528                 update_caldac(cdac);
 529                 //usleep(100000);
 530
 531                 new_sv_measure(&sv);
 532
 533                 l->y_data[i]=sv.average;
 534                 if(!isnan(sv.average)){
 535                         sum_err+=sv.error;
 536                         sum_err_count++;
 537                 }
 538                 l->n++;
 539         }
 540
 541         caldacs[cdac].current=orig;
 542         update_caldac(cdac);
 543
 544         l->yerr=sum_err/sum_err_count;
 545         l->dx=step;
 546         l->x0=0;
 547
 548         linear_fit_monotonic(l);
 549
 550         if(verbose>=2 || (verbose>=1 && fabs(l->slope/l->err_slope)>4.0)){
 551                 sci_sprint_alt(str,l->slope,l->err_slope);
 552                 printf("caldac[%d] gain=%s V/bit S_min=%g dof=%g\n",
 553                         cdac,str,l->S_min,l->dof);
 554                 //printf("--> %g\n",fabs(l.slope/l.err_slope));
 555         }
 556
 557         if(verbose>=3)dump_curve(l);
 558
 559         free(l->y_data);
 560
 561         return l->slope;
 562 }
 563
 564
 565 double check_gain_chan_fine(linear_fit_t *l,unsigned int ad_chanspec,int cdac)
 566 {
 567         int orig,i,n;
 568         int step;
 569         new_sv_t sv;
 570         double sum_err;
 571         int sum_err_count=0;
 572         char str[20];
 573         int fine_size = 10;
 574
 575         n=2*fine_size+1;
 576         memset(l,0,sizeof(*l));
 577
 578         step=1;
 579         l->n=0;
 580
 581         l->y_data=malloc(n*sizeof(double)/step);
 582         if(l->y_data == NULL)
 583         {
 584                 perror("comedi_calibrate");
 585                 exit(1);
 586         }
 587
 588         orig=caldacs[cdac].current;
 589
 590         new_sv_init(&sv,dev,0,
 591                 CR_CHAN(ad_chanspec),
 592                 CR_RANGE(ad_chanspec),
 593                 CR_AREF(ad_chanspec));
 594
 595         caldacs[cdac].current=0;
 596         update_caldac(cdac);
 597         usleep(100000);
 598
 599         new_sv_measure(&sv);
 600
 601         sum_err=0;
 602         for(i=0;i<n;i++){
 603                 caldacs[cdac].current=i+orig-fine_size;
 604                 update_caldac(cdac);
 605                 usleep(100000);
 606
 607                 new_sv_measure(&sv);
 608
 609                 l->y_data[i]=sv.average;
 610                 if(!isnan(sv.average)){
 611                         sum_err+=sv.error;
 612                         sum_err_count++;
 613                 }
 614                 l->n++;
 615         }
 616
 617         caldacs[cdac].current=orig;
 618         update_caldac(cdac);
 619
 620         l->yerr=sum_err/sum_err_count;
 621         l->dx=1;
 622         l->x0=orig-fine_size;
 623
 624         linear_fit_monotonic(l);
 625
 626         if(verbose>=2 || (verbose>=1 && fabs(l->slope/l->err_slope)>4.0)){
 627                 sci_sprint_alt(str,l->slope,l->err_slope);
 628                 printf("caldac[%d] gain=%s V/bit S_min=%g dof=%g\n",
 629                         cdac,str,l->S_min,l->dof);
 630                 //printf("--> %g\n",fabs(l.slope/l.err_slope));
 631         }
 632
 633         if(verbose>=3)dump_curve(l);
 634
 635         free(l->y_data);
 636
 637         return l->slope;
 638 }
 639
 640
 641
 642 /* helpers */
 643
 644 int get_bipolar_lowgain(comedi_t *dev,int subdev)
 645 {
 646         int ret = -1;
 647         int i;
 648         int n_ranges = comedi_get_n_ranges(dev,subdev,0);
 649         double max = 0;
 650         comedi_range *range;
 651
 652         for(i=0;i<n_ranges;i++){
 653                 range = comedi_get_range(dev,subdev,0,i);
 654                 if(range->min != -range->max)continue;
 655                 if(range->max>max){
 656                         ret = i;
 657                         max=range->max;
 658                 }
 659         }
 660
 661         return ret;
 662 }
 663
 664 int get_bipolar_highgain(comedi_t *dev,int subdev)
 665 {
 666         int ret = -1;
 667         int i;
 668         int n_ranges = comedi_get_n_ranges(dev,subdev,0);
 669         double min = HUGE_VAL;
 670         comedi_range *range;
 671
 672         for(i=0;i<n_ranges;i++){
 673                 range = comedi_get_range(dev,subdev,0,i);
 674                 if(range->min != -range->max)continue;
 675                 if(range->max<min){
 676                         ret = i;
 677                         min=range->max;
 678                 }
 679         }
 680
 681         return ret;
 682 }
 683
 684 int get_unipolar_lowgain(comedi_t *dev,int subdev)
 685 {
 686         int ret = -1;
 687         int i;
 688         int n_ranges = comedi_get_n_ranges(dev,subdev,0);
 689         double max = 0;
 690         comedi_range *range;
 691
 692         for(i=0;i<n_ranges;i++){
 693                 range = comedi_get_range(dev,subdev,0,i);
 694                 if(range->min != 0)continue;
 695                 if(range->max>max){
 696                         ret = i;
 697                         max=range->max;
 698                 }
 699         }
 700
 701         return ret;
 702 }
 703
 704
 705 int read_eeprom(int addr)
 706 {
 707         unsigned int data=0;
 708
 709         comedi_data_read(dev,eeprom_subdev,addr,0,0,&data);
 710
 711         return data;
 712 }
 713
 714 double read_chan(int adc,int range)
 715 {
 716         int n;
 717         new_sv_t sv;
 718         char str[20];
 719
 720         new_sv_init(&sv,dev,0,adc,range,AREF_OTHER);
 721         sv.order=7;
 722         n=new_sv_measure(&sv);
 723
 724         sci_sprint_alt(str,sv.average,sv.error);
 725         printf("chan=%d ave=%s\n",adc,str);
 726
 727         return sv.average;
 728 }
 729
 730 int read_chan2(char *s,int adc,int range)
 731 {
 732         int n;
 733         new_sv_t sv;
 734
 735         new_sv_init(&sv,dev,0,adc,range,AREF_OTHER);
 736         sv.order=7;
 737         n=new_sv_measure(&sv);
 738
 739         return sci_sprint_alt(s,sv.average,sv.error);
 740 }
 741
 742 #if 0
 743 void set_ao(comedi_t *dev,int subdev,int chan,int range,double value)
 744 {
 745         comedi_range *r = comedi_get_range(dev,subdev,chan,range);
 746         lsampl_t maxdata = comedi_get_maxdata(dev,subdev,chan);
 747         lsampl_t data;
 748
 749         data = comedi_from_phys(value,r,maxdata);
 750
 751         comedi_data_write(dev,subdev,chan,range,AREF_GROUND,data);
 752 }
 753 #endif
 754
 755
 756 int new_sv_init(new_sv_t *sv,comedi_t *dev,int subdev,int chan,int range,int aref)
 757 {
 758         memset(sv,0,sizeof(*sv));
 759
 760         sv->subd=subdev;
 761         //sv->t.flags=TRIG_DITHER;
 762         sv->chan=chan;
 763         sv->range=range;
 764         sv->aref=aref;
 765
 766         //sv->chanlist[0]=CR_PACK(chan,range,aref);
 767
 768         sv->maxdata=comedi_get_maxdata(dev,subdev,chan);
 769         sv->rng=comedi_get_range(dev,subdev,chan,range);
 770
 771         sv->order=7;
 772
 773         return 0;
 774 }
 775
 776 int comedi_data_read_n(comedi_t *it,unsigned int subdev,unsigned int chan,
 777         unsigned int range,unsigned int aref,lsampl_t *data,unsigned int n)
 778 {
 779         comedi_insn insn;
 780         int ret;
 781
 782         if(n==0)return 0;
 783
 784         insn.insn = INSN_READ;
 785         insn.n = n;
 786         insn.data = data;
 787         insn.subdev = subdev;
 788         insn.chanspec = CR_PACK(chan,range,aref);
 789         /* enable dithering */
 790         insn.chanspec |= (1<<26);
 791
 792         ret = comedi_do_insn(it,&insn);
 793
 794         if(ret>0)return n;
 795
 796         printf("insn barfed: subdev=%d, chan=%d, range=%d, aref=%d, "
 797                 "n=%d, ret=%d, %s\n",subdev,chan,range,aref,n,ret,
 798                 strerror(errno));
 799
 800         return ret;
 801 }
 802
 803 int new_sv_measure(new_sv_t *sv)
 804 {
 805         lsampl_t *data;
 806         int n,i,ret;
 807
 808         double a,x,s,s2;
 809
 810         n=1<<sv->order;
 811
 812         data=malloc(sizeof(lsampl_t)*n);
 813         if(data == NULL)
 814         {
 815                 perror("comedi_calibrate");
 816                 exit(1);
 817         }
 818
 819         for(i=0;i<n;){
 820                 ret = comedi_data_read_n(dev,sv->subd,sv->chan,sv->range,
 821                         sv->aref,data+i,n-i);
 822                 if(ret<0){
 823                         printf("barf\n");
 824                         goto out;
 825                 }
 826                 i+=ret;
 827         }
 828
 829         s=0;
 830         s2=0;
 831         a=comedi_to_phys(data[0],sv->rng,sv->maxdata);
 832         for(i=0;i<n;i++){
 833                 x=comedi_to_phys(data[i],sv->rng,sv->maxdata);
 834                 s+=x-a;
 835                 s2+=(x-a)*(x-a);
 836         }
 837         sv->average=a+s/n;
 838         sv->stddev=sqrt(n*s2-s*s)/n;
 839         sv->error=sv->stddev/sqrt(n);
 840
 841         ret=n;
 842
 843 out:
 844         free(data);
 845
 846         return ret;
 847 }
 848
 849 int new_sv_measure_order(new_sv_t *sv,int order)
 850 {
 851         lsampl_t *data;
 852         int n,i,ret;
 853         double a,x,s,s2;
 854
 855         n=1<<order;
 856
 857         data=malloc(sizeof(lsampl_t)*n);
 858         if(data == NULL)
 859         {
 860                 perror("comedi_calibrate");
 861                 exit(1);
 862         }
 863
 864         for(i=0;i<n;){
 865                 ret = comedi_data_read_n(dev,sv->subd,sv->chan,sv->range,
 866                         sv->aref,data+i,n-i);
 867                 if(ret<0){
 868                         printf("barf order\n");
 869                         goto out;
 870                 }
 871                 i+=ret;
 872         }
 873
 874         s=0;
 875         s2=0;
 876         a=comedi_to_phys(data[0],sv->rng,sv->maxdata);
 877         for(i=0;i<n;i++){
 878                 x=comedi_to_phys(data[i],sv->rng,sv->maxdata);
 879                 s+=x-a;
 880                 s2+=(x-a)*(x-a);
 881         }
 882         sv->average=a+s/n;
 883         sv->stddev=sqrt(n*s2-s*s)/n;
 884         sv->error=sv->stddev/sqrt(n);
 885
 886         ret=n;
 887
 888 out:
 889         free(data);
 890
 891         return ret;
 892 }
 893
 894
 895
 896 /* linear fitting */
 897
 898 int calculate_residuals(linear_fit_t *l);
 899
 900 int linear_fit_monotonic(linear_fit_t *l)
 901 {
 902         double x,y;
 903         double sxp;
 904         int i;
 905
 906         l->min=HUGE_VAL;
 907         l->max=-HUGE_VAL;
 908         l->s1=0;
 909         l->sx=0;
 910         l->sy=0;
 911         l->sxy=0;
 912         l->sxx=0;
 913         for(i=0;i<l->n;i++){
 914                 x=l->x0+i*l->dx;
 915                 y=l->y_data[i];
 916
 917                 if(isnan(y))continue;
 918
 919                 if(l->y_data[i]<l->min)l->min=l->y_data[i];
 920                 if(l->y_data[i]>l->max)l->max=l->y_data[i];
 921                 l->s1+=1;
 922                 l->sx+=x;
 923                 l->sy+=y;
 924                 l->sxy+=x*y;
 925                 l->sxx+=x*x;
 926         }
 927         sxp=l->sxx-l->sx*l->sx/l->s1;
 928
 929         l->ave_x=l->sx/l->s1;
 930         l->ave_y=l->sy/l->s1;
 931         l->slope=(l->s1*l->sxy-l->sx*l->sy)/(l->s1*l->sxx-l->sx*l->sx);
 932         l->err_slope=l->yerr/sqrt(sxp);
 933         l->err_ave_y=l->yerr/sqrt(l->s1);
 934
 935         calculate_residuals(l);
 936
 937         return 0;
 938 }
 939
 940 int calculate_residuals(linear_fit_t *l)
 941 {
 942         double x,y;
 943         double res,sum_res2;
 944         int i;
 945
 946         sum_res2=0;
 947         for(i=0;i<l->n;i++){
 948                 x=l->x0+i*l->dx-l->ave_x;
 949                 y=l->y_data[i];
 950
 951                 if(isnan(y))continue;
 952
 953                 res=l->ave_y+l->slope*x-y;
 954
 955                 sum_res2+=res*res;
 956         }
 957         l->S_min=sum_res2/(l->yerr*l->yerr);
 958         l->dof=l->s1-2;
 959
 960         return 0;
 961 }
 962
 963 double linear_fit_func_y(linear_fit_t *l,double x)
 964 {
 965         return l->ave_y+l->slope*(x-l->ave_x);
 966 }
 967
 968 double linear_fit_func_x(linear_fit_t *l,double y)
 969 {
 970         return l->ave_x+(y-l->ave_y)/l->slope;
 971 }
 972
 973 void dump_curve(linear_fit_t *l)
 974 {
 975         static int dump_number=0;
 976         double x,y;
 977         int i;
 978
 979         printf("start dump %d\n",dump_number);
 980         for(i=0;i<l->n;i++){
 981                 x=l->x0+i*l->dx-l->ave_x;
 982                 y=l->y_data[i];
 983                 printf("D%d: %d %g %g %g\n",dump_number,i,y,
 984                         l->ave_y+l->slope*x,
 985                         l->ave_y+l->slope*x-y);
 986         }
 987         printf("end dump\n");
 988         dump_number++;
 989 }
 990
 991
 992 /* printing of scientific numbers (with errors) */
 993
 994 int sci_sprint(char *s,double x,double y)
 995 {
 996         int errsig;
 997         int maxsig;
 998         int sigfigs;
 999         double mantissa;
1000         double error;
1001         double mindigit;
1002
1003         errsig = floor(log10(y));
1004         maxsig = floor(log10(x));
1005         mindigit = pow(10,errsig);
1006
1007         if(maxsig<errsig)maxsig=errsig;
1008
1009         sigfigs = maxsig-errsig+2;
1010
1011         mantissa = x*pow(10,-maxsig);
1012         error = y*pow(10,-errsig+1);
1013
1014         return sprintf(s,"%0.*f(%2.0f)e%d",sigfigs-1,mantissa,error,maxsig);
1015 }
1016
1017 int sci_sprint_alt(char *s,double x,double y)
1018 {
1019         int errsig;
1020         int maxsig;
1021         int sigfigs;
1022         double mantissa;
1023         double error;
1024         double mindigit;
1025
1026         errsig = floor(log10(fabs(y)));
1027         maxsig = floor(log10(fabs(x)));
1028         mindigit = pow(10,errsig);
1029
1030         if(maxsig<errsig)maxsig=errsig;
1031
1032         sigfigs = maxsig-errsig+2;
1033
1034         mantissa = x*pow(10,-maxsig);
1035         error = y*pow(10,-errsig+1);
1036
1037         if(isnan(x)){
1038                 return sprintf(s,"%g",x);
1039         }
1040         if(errsig==1 && maxsig<4 && maxsig>1){
1041                 return sprintf(s,"%0.0f(%2.0f)",x,error);
1042         }
1043         if(maxsig<=0 && maxsig>=-2){
1044                 return sprintf(s,"%0.*f(%2.0f)",sigfigs-1-maxsig,
1045                         mantissa*pow(10,maxsig),error);
1046         }
1047         return sprintf(s,"%0.*f(%2.0f)e%d",sigfigs-1,mantissa,error,maxsig);
1048 }
1049