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 'f':
 114                         fn = optarg;
 115                         break;
 116                 case 'v':
 117                         verbose++;
 118                         break;
 119                 case 'q':
 120                         verbose--;
 121                         break;
 122                 case 0x1000:
 123                         drivername = optarg;
 124                         break;
 125                 case 0x1001:
 126                         devicename = optarg;
 127                         break;
 128                 default:
 129                         printf("bad option\n");
 130                         exit(1);
 131                 }
 132         }
 133
 134         dev = comedi_open(fn);
 135         if (dev == NULL ) {
 136                 perror(fn);
 137                 exit(0);
 138         }
 139
 140         if(!drivername)
 141                 drivername=comedi_get_driver_name(dev);
 142         if(!devicename)
 143                 devicename=comedi_get_board_name(dev);
 144
 145         ad_subdev=comedi_find_subdevice_by_type(dev,COMEDI_SUBD_AI,0);
 146         da_subdev=comedi_find_subdevice_by_type(dev,COMEDI_SUBD_AO,0);
 147         caldac_subdev=comedi_find_subdevice_by_type(dev,COMEDI_SUBD_CALIB,0);
 148         eeprom_subdev=comedi_find_subdevice_by_type(dev,COMEDI_SUBD_MEMORY,0);
 149
 150         for(i=0;i<n_drivers;i++){
 151                 if(!strcmp(drivers[i].name,drivername)){
 152                         this_board = drivers+i;
 153                         goto ok;
 154                 }
 155         }
 156         printf("Driver %s unknown\n",drivername);
 157         return 1;
 158
 159 ok:
 160         this_board->setup();
 161
 162         if(device_status<STATUS_DONE){
 163                 printf("Warning: device not fully calibrated due to insufficient information\n");
 164                 printf("Please send this output to <ds@schleef.org>\n");
 165                 if(verbose<0)verbose=0;
 166                 if(device_status==STATUS_UNKNOWN){
 167                         do_reset=1;
 168                         do_dump=1;
 169                         do_calibrate=0;
 170                         do_results=0;
 171                 }
 172                 if(device_status==STATUS_SOME){
 173                         do_reset=1;
 174                         do_dump=1;
 175                         do_calibrate=1;
 176                         do_results=1;
 177                 }
 178         }
 179         if(verbose>=0){
 180                 printf("$Id$\n");
 181                 printf("Driver name: %s\n",drivername);
 182                 printf("Device name: %s\n",devicename);
 183                 printf("Comedi version: %d.%d.%d\n",
 184                         (comedi_get_version_code(dev)>>16)&0xff,
 185                         (comedi_get_version_code(dev)>>8)&0xff,
 186                         (comedi_get_version_code(dev))&0xff);
 187         }
 188
 189         if(do_reset)reset_caldacs();
 190         if(do_dump)observe();
 191         if(do_calibrate && do_cal)do_cal();
 192         if(do_results)observe();
 193
 194         return 0;
 195 }
 196
 197 void set_target(int obs,double target)
 198 {
 199         comedi_range *range;
 200         lsampl_t maxdata, data;
 201
 202         range = comedi_get_range(dev,observables[obs].preobserve_insn.subdev,
 203                 CR_CHAN(observables[obs].preobserve_insn.chanspec),
 204                 CR_RANGE(observables[obs].preobserve_insn.chanspec));
 205         maxdata = comedi_get_maxdata(dev,observables[obs].preobserve_insn.subdev,
 206                 CR_CHAN(observables[obs].preobserve_insn.chanspec));
 207
 208         data = comedi_from_phys(target,range,maxdata);
 209
 210         observables[obs].preobserve_data = data;
 211         observables[obs].target = comedi_to_phys(data,range,maxdata);
 212 }
 213
 214 void observe(void)
 215 {
 216         int i;
 217
 218         for(i=0;i<n_observables;i++){
 219                 preobserve(i);
 220                 DPRINT(0,"%s\n",observables[i].name);
 221                 measure_observable(i);
 222                 observable_dependence(i);
 223         }
 224
 225 }
 226
 227 void preobserve(int obs)
 228 {
 229         if(observables[obs].preobserve_insn.n!=0){
 230                 comedi_do_insn(dev,&observables[obs].preobserve_insn);
 231         }
 232 }
 233
 234 void measure_observable(int obs)
 235 {
 236         char s[32];
 237         int n;
 238         new_sv_t sv;
 239
 240         new_sv_init(&sv,dev,
 241                 observables[obs].observe_insn.subdev,
 242                 CR_CHAN(observables[obs].observe_insn.chanspec),
 243                 CR_RANGE(observables[obs].observe_insn.chanspec),
 244                 CR_AREF(observables[obs].observe_insn.chanspec));
 245         sv.order=7;
 246         n=new_sv_measure(&sv);
 247
 248         sci_sprint_alt(s,sv.average,sv.error);
 249         printf("offset %s, target %g\n",s,observables[obs].target);
 250 }
 251
 252 void observable_dependence(int obs)
 253 {
 254         int i;
 255         linear_fit_t l;
 256
 257         for(i=0;i<n_caldacs;i++){
 258                 check_gain_chan_x(&l,
 259                         observables[obs].observe_insn.chanspec,i);
 260         }
 261
 262 }
 263
 264
 265 void postgain_cal(int obs1, int obs2, int dac)
 266 {
 267         double offset1,offset2;
 268         linear_fit_t l;
 269         double slope1,slope2;
 270         double a;
 271         double gain;
 272         comedi_range *range1,*range2;
 273
 274         DPRINT(0,"postgain calibration\n");
 275         preobserve(obs1);
 276         check_gain_chan_x(&l,observables[obs1].observe_insn.chanspec,dac);
 277         offset1=linear_fit_func_y(&l,caldacs[dac].current);
 278         DPRINT(1,"obs1: [%d] offset %g\n",obs1,offset1);
 279         range1 = comedi_get_range(dev,observables[obs1].observe_insn.subdev,
 280                 CR_CHAN(observables[obs1].observe_insn.chanspec),
 281                 CR_RANGE(observables[obs1].observe_insn.chanspec));
 282         slope1=l.slope;
 283
 284         preobserve(obs2);
 285         check_gain_chan_x(&l,observables[obs2].observe_insn.chanspec,dac);
 286         offset2=linear_fit_func_y(&l,caldacs[dac].current);
 287         DPRINT(1,"obs2: [%d] offset %g\n",obs2,offset2);
 288         range2 = comedi_get_range(dev,observables[obs2].observe_insn.subdev,
 289                 CR_CHAN(observables[obs2].observe_insn.chanspec),
 290                 CR_RANGE(observables[obs2].observe_insn.chanspec));
 291         slope2=l.slope;
 292
 293         gain = (range1->max-range1->min)/(range2->max-range2->min);
 294         DPRINT(3,"range1 %g range2 %g\n", range1->max-range1->min,
 295                 range2->max-range2->min);
 296         DPRINT(2,"gain: %g\n",gain);
 297
 298         DPRINT(2,"difference: %g\n",offset2-offset1);
 299
 300         a = (offset1-offset2)/(slope1-slope2);
 301         a=caldacs[dac].current-a;
 302
 303         DPRINT(0,"caldac[%d] set to %g\n",dac,a);
 304
 305         caldacs[dac].current=rint(a);
 306         update_caldac(dac);
 307         usleep(100000);
 308
 309         if(verbose>=2){
 310                 preobserve(obs1);
 311                 measure_observable(obs1);
 312                 preobserve(obs2);
 313                 measure_observable(obs2);
 314         }
 315 }
 316
 317 void cal1(int obs, int dac)
 318 {
 319         linear_fit_t l;
 320         double offset;
 321         double target;
 322         double gain;
 323         double a;
 324
 325         DPRINT(0,"cal1\n");
 326         preobserve(obs);
 327         check_gain_chan_x(&l,observables[obs].observe_insn.chanspec,dac);
 328         offset=linear_fit_func_y(&l,caldacs[dac].current);
 329         gain=l.slope;
 330
 331         target = observables[obs].target;
 332         a=caldacs[dac].current+(target-offset)/gain;
 333
 334         caldacs[dac].current=rint(a);
 335         update_caldac(dac);
 336         usleep(100000);
 337
 338         DPRINT(0,"caldac[%d] set to %g\n",dac,a);
 339         if(verbose>=2){
 340                 measure_observable(obs);
 341         }
 342 }
 343
 344 void chan_cal(int adc,int cdac,int range,double target)
 345 {
 346         linear_fit_t l;
 347         double offset;
 348         double gain;
 349         double a;
 350         char s[32];
 351
 352         check_gain_chan_x(&l,CR_PACK(adc,range,AREF_OTHER),cdac);
 353         offset=linear_fit_func_y(&l,caldacs[cdac].current);
 354         gain=l.slope;
 355
 356         a=caldacs[cdac].current+(target-offset)/gain;
 357
 358         caldacs[cdac].current=rint(a);
 359         update_caldac(cdac);
 360
 361         read_chan2(s,adc,range);
 362         DPRINT(0,"caldac[%d] set to %g, offset=%s\n",cdac,a,s);
 363 }
 364
 365
 366
 367 void channel_dependence(int adc,int range)
 368 {
 369         int i;
 370         double gain;
 371
 372         for(i=0;i<n_caldacs;i++){
 373                 gain=check_gain_chan(adc,range,i);
 374         }
 375 }
 376
 377 void caldac_dependence(int caldac)
 378 {
 379         int i;
 380         double gain;
 381
 382         for(i=0;i<8;i++){
 383                 gain=check_gain_chan(i,0,caldac);
 384         }
 385 }
 386
 387 void dump_curve(int adc,int caldac)
 388 {
 389         linear_fit_t l;
 390
 391         check_gain_chan_x(&l,CR_PACK(adc,0,AREF_OTHER),caldac);
 392 }
 393
 394
 395 void setup_caldacs(void)
 396 {
 397         int n_chan,i;
 398
 399         if(caldac_subdev<0){
 400                 printf("no calibration subdevice\n");
 401                 return;
 402         }
 403
 404         n_chan=comedi_get_n_channels(dev,caldac_subdev);
 405
 406         for(i=0;i<n_chan;i++){
 407                 caldacs[i].subdev=caldac_subdev;
 408                 caldacs[i].chan=i;
 409                 caldacs[i].maxdata=comedi_get_maxdata(dev,caldac_subdev,i);
 410                 caldacs[i].current=0;
 411         }
 412
 413         n_caldacs=n_chan;
 414 }
 415
 416 void reset_caldacs(void)
 417 {
 418         int i;
 419
 420         for(i=0;i<n_caldacs;i++){
 421                 caldacs[i].current=caldacs[i].maxdata/2;
 422                 update_caldac(i);
 423         }
 424 }
 425
 426 void update_caldac(int i)
 427 {
 428         int ret;
 429
 430         DPRINT(3,"update %d %d %d\n",caldacs[i].subdev,caldacs[i].chan,caldacs[i].current);
 431         if(caldacs[i].current<0){
 432                 DPRINT(0,"caldac set out of range (%d<0)\n",caldacs[i].current);
 433                 caldacs[i].current=0;
 434         }
 435         if(caldacs[i].current>caldacs[i].maxdata){
 436                 DPRINT(0,"caldac set out of range (%d>%d)\n",
 437                         caldacs[i].current,caldacs[i].maxdata);
 438                 caldacs[i].current=caldacs[i].maxdata;
 439         }
 440
 441         ret = comedi_data_write(dev,caldacs[i].subdev,caldacs[i].chan,0,0,
 442                 caldacs[i].current);
 443         if(ret<0)perror("update_caldac()");
 444 }
 445
 446
 447 void check_gain(int ad_chan,int range)
 448 {
 449         int i;
 450
 451         for(i=0;i<n_caldacs;i++){
 452                 check_gain_chan(ad_chan,range,i);
 453         }
 454 }
 455
 456 double check_gain_chan(int ad_chan,int range,int cdac)
 457 {
 458         linear_fit_t l;
 459
 460         return check_gain_chan_x(&l,CR_PACK(ad_chan,range,AREF_OTHER),cdac);
 461 }
 462
 463 double check_gain_chan_x(linear_fit_t *l,unsigned int ad_chanspec,int cdac)
 464 {
 465         int orig,i,n;
 466         int step;
 467         new_sv_t sv;
 468         double sum_err;
 469         int sum_err_count=0;
 470         char str[20];
 471
 472         n=caldacs[cdac].maxdata+1;
 473         memset(l,0,sizeof(*l));
 474
 475         step=n/256;
 476         if(step<1)step=1;
 477         l->n=0;
 478
 479         l->y_data=malloc(n*sizeof(double)/step);
 480         if(l->y_data == NULL)
 481         {
 482                 perror("comedi_calibrate");
 483                 exit(1);
 484         }
 485
 486         orig=caldacs[cdac].current;
 487
 488         new_sv_init(&sv,dev,0,
 489                 CR_CHAN(ad_chanspec),
 490                 CR_RANGE(ad_chanspec),
 491                 CR_AREF(ad_chanspec));
 492
 493         caldacs[cdac].current=0;
 494         update_caldac(cdac);
 495
 496         new_sv_measure(&sv);
 497         usleep(100000);
 498
 499         sum_err=0;
 500         for(i=0;i*step<n;i++){
 501                 caldacs[cdac].current=i*step;
 502                 update_caldac(cdac);
 503
 504                 new_sv_measure(&sv);
 505
 506                 l->y_data[i]=sv.average;
 507                 if(!isnan(sv.average)){
 508                         sum_err+=sv.error;
 509                         sum_err_count++;
 510                 }
 511                 l->n++;
 512         }
 513
 514         caldacs[cdac].current=orig;
 515         update_caldac(cdac);
 516
 517         l->yerr=sum_err/sum_err_count;
 518         l->dx=step;
 519         l->x0=0;
 520
 521         linear_fit_monotonic(l);
 522
 523         if(verbose>=1 || (verbose>=0 && fabs(l->slope/l->err_slope)>4.0)){
 524                 sci_sprint_alt(str,l->slope,l->err_slope);
 525                 printf("caldac[%d] gain=%s V/bit S_min=%g dof=%g\n",
 526                         cdac,str,l->S_min,l->dof);
 527                 //printf("--> %g\n",fabs(l.slope/l.err_slope));
 528         }
 529
 530         if(verbose>=2){
 531                 static int dump_number=0;
 532                 double x,y;
 533
 534                 printf("start dump %d\n",dump_number);
 535                 for(i=0;i<l->n;i++){
 536                         x=l->x0+i*l->dx-l->ave_x;
 537                         y=l->y_data[i];
 538                         printf("D%d: %d %g %g %g\n",dump_number,i,y,
 539                                 l->ave_y+l->slope*x,
 540                                 l->ave_y+l->slope*x-y);
 541                 }
 542                 printf("end dump\n");
 543                 dump_number++;
 544         }
 545
 546
 547         free(l->y_data);
 548
 549         return l->slope;
 550 }
 551
 552
 553
 554
 555 /* helpers */
 556
 557 int get_bipolar_lowgain(comedi_t *dev,int subdev)
 558 {
 559         int ret = -1;
 560         int i;
 561         int n_ranges = comedi_get_n_ranges(dev,subdev,0);
 562         double max = 0;
 563         comedi_range *range;
 564
 565         for(i=0;i<n_ranges;i++){
 566                 range = comedi_get_range(dev,subdev,0,i);
 567                 if(range->min != -range->max)continue;
 568                 if(range->max>max){
 569                         ret = i;
 570                         max=range->max;
 571                 }
 572         }
 573
 574         return ret;
 575 }
 576
 577 int get_bipolar_highgain(comedi_t *dev,int subdev)
 578 {
 579         int ret = -1;
 580         int i;
 581         int n_ranges = comedi_get_n_ranges(dev,subdev,0);
 582         double min = HUGE_VAL;
 583         comedi_range *range;
 584
 585         for(i=0;i<n_ranges;i++){
 586                 range = comedi_get_range(dev,subdev,0,i);
 587                 if(range->min != -range->max)continue;
 588                 if(range->max<min){
 589                         ret = i;
 590                         min=range->max;
 591                 }
 592         }
 593
 594         return ret;
 595 }
 596
 597 int get_unipolar_lowgain(comedi_t *dev,int subdev)
 598 {
 599         int ret = -1;
 600         int i;
 601         int n_ranges = comedi_get_n_ranges(dev,subdev,0);
 602         double max = 0;
 603         comedi_range *range;
 604
 605         for(i=0;i<n_ranges;i++){
 606                 range = comedi_get_range(dev,subdev,0,i);
 607                 if(range->min != 0)continue;
 608                 if(range->max>max){
 609                         ret = i;
 610                         max=range->max;
 611                 }
 612         }
 613
 614         return ret;
 615 }
 616
 617
 618 int read_eeprom(int addr)
 619 {
 620         unsigned int data=0;
 621
 622         comedi_data_read(dev,eeprom_subdev,addr,0,0,&data);
 623
 624         return data;
 625 }
 626
 627 double read_chan(int adc,int range)
 628 {
 629         int n;
 630         new_sv_t sv;
 631         char str[20];
 632
 633         new_sv_init(&sv,dev,0,adc,range,AREF_OTHER);
 634         sv.order=7;
 635         n=new_sv_measure(&sv);
 636
 637         sci_sprint_alt(str,sv.average,sv.error);
 638         printf("chan=%d ave=%s\n",adc,str);
 639
 640         return sv.average;
 641 }
 642
 643 int read_chan2(char *s,int adc,int range)
 644 {
 645         int n;
 646         new_sv_t sv;
 647
 648         new_sv_init(&sv,dev,0,adc,range,AREF_OTHER);
 649         sv.order=7;
 650         n=new_sv_measure(&sv);
 651
 652         return sci_sprint_alt(s,sv.average,sv.error);
 653 }
 654
 655 void set_ao(comedi_t *dev,int subdev,int chan,int range,double value)
 656 {
 657         comedi_range *r = comedi_get_range(dev,subdev,chan,range);
 658         lsampl_t maxdata = comedi_get_maxdata(dev,subdev,chan);
 659         lsampl_t data;
 660
 661         data = comedi_from_phys(value,r,maxdata);
 662
 663         comedi_data_write(dev,subdev,chan,range,AREF_GROUND,data);
 664 }
 665
 666
 667 int new_sv_init(new_sv_t *sv,comedi_t *dev,int subdev,int chan,int range,int aref)
 668 {
 669         memset(sv,0,sizeof(*sv));
 670
 671         sv->subd=subdev;
 672         //sv->t.flags=TRIG_DITHER;
 673         sv->chan=chan;
 674         sv->range=range;
 675         sv->aref=aref;
 676
 677         //sv->chanlist[0]=CR_PACK(chan,range,aref);
 678
 679         sv->maxdata=comedi_get_maxdata(dev,subdev,chan);
 680         sv->rng=comedi_get_range(dev,subdev,chan,range);
 681
 682         sv->order=7;
 683
 684         return 0;
 685 }
 686
 687 int comedi_data_read_n(comedi_t *it,unsigned int subdev,unsigned int chan,
 688         unsigned int range,unsigned int aref,lsampl_t *data,unsigned int n)
 689 {
 690         comedi_insn insn;
 691         int ret;
 692
 693         if(n==0)return 0;
 694
 695         insn.insn = INSN_READ;
 696         insn.n = n;
 697         insn.data = data;
 698         insn.subdev = subdev;
 699         insn.chanspec = CR_PACK(chan,range,aref);
 700         /* enable dithering */
 701         insn.chanspec |= (1<<26);
 702
 703         ret = comedi_do_insn(it,&insn);
 704
 705         if(ret>0)return n;
 706
 707         printf("insn barfed: subdev=%d, chan=%d, range=%d, aref=%d, "
 708                 "n=%d, ret=%d, %s\n",subdev,chan,range,aref,n,ret,
 709                 strerror(errno));
 710
 711         return ret;
 712 }
 713
 714 int new_sv_measure(new_sv_t *sv)
 715 {
 716         lsampl_t *data;
 717         int n,i,ret;
 718
 719         double a,x,s,s2;
 720
 721         n=1<<sv->order;
 722
 723         data=malloc(sizeof(lsampl_t)*n);
 724         if(data == NULL)
 725         {
 726                 perror("comedi_calibrate");
 727                 exit(1);
 728         }
 729
 730         for(i=0;i<n;){
 731                 ret = comedi_data_read_n(dev,sv->subd,sv->chan,sv->range,
 732                         sv->aref,data+i,n-i);
 733                 if(ret<0){
 734                         printf("barf\n");
 735                         goto out;
 736                 }
 737                 i+=ret;
 738         }
 739
 740         s=0;
 741         s2=0;
 742         a=comedi_to_phys(data[0],sv->rng,sv->maxdata);
 743         for(i=0;i<n;i++){
 744                 x=comedi_to_phys(data[i],sv->rng,sv->maxdata);
 745                 s+=x-a;
 746                 s2+=(x-a)*(x-a);
 747         }
 748         sv->average=a+s/n;
 749         sv->stddev=sqrt(n*s2-s*s)/n;
 750         sv->error=sv->stddev/sqrt(n);
 751
 752         ret=n;
 753
 754 out:
 755         free(data);
 756
 757         return ret;
 758 }
 759
 760 int new_sv_measure_order(new_sv_t *sv,int order)
 761 {
 762         lsampl_t *data;
 763         int n,i,ret;
 764         double a,x,s,s2;
 765
 766         n=1<<order;
 767
 768         data=malloc(sizeof(lsampl_t)*n);
 769         if(data == NULL)
 770         {
 771                 perror("comedi_calibrate");
 772                 exit(1);
 773         }
 774
 775         for(i=0;i<n;){
 776                 ret = comedi_data_read_n(dev,sv->subd,sv->chan,sv->range,
 777                         sv->aref,data+i,n-i);
 778                 if(ret<0){
 779                         printf("barf order\n");
 780                         goto out;
 781                 }
 782                 i+=ret;
 783         }
 784
 785         s=0;
 786         s2=0;
 787         a=comedi_to_phys(data[0],sv->rng,sv->maxdata);
 788         for(i=0;i<n;i++){
 789                 x=comedi_to_phys(data[i],sv->rng,sv->maxdata);
 790                 s+=x-a;
 791                 s2+=(x-a)*(x-a);
 792         }
 793         sv->average=a+s/n;
 794         sv->stddev=sqrt(n*s2-s*s)/n;
 795         sv->error=sv->stddev/sqrt(n);
 796
 797         ret=n;
 798
 799 out:
 800         free(data);
 801
 802         return ret;
 803 }
 804
 805
 806
 807 /* linear fitting */
 808
 809 int calculate_residuals(linear_fit_t *l);
 810
 811 int linear_fit_monotonic(linear_fit_t *l)
 812 {
 813         double x,y;
 814         double sxp;
 815         int i;
 816
 817         l->min=HUGE_VAL;
 818         l->max=-HUGE_VAL;
 819         l->s1=0;
 820         l->sx=0;
 821         l->sy=0;
 822         l->sxy=0;
 823         l->sxx=0;
 824         for(i=0;i<l->n;i++){
 825                 x=l->x0+i*l->dx;
 826                 y=l->y_data[i];
 827
 828                 if(isnan(y))continue;
 829
 830                 if(l->y_data[i]<l->min)l->min=l->y_data[i];
 831                 if(l->y_data[i]>l->max)l->max=l->y_data[i];
 832                 l->s1+=1;
 833                 l->sx+=x;
 834                 l->sy+=y;
 835                 l->sxy+=x*y;
 836                 l->sxx+=x*x;
 837         }
 838         sxp=l->sxx-l->sx*l->sx/l->s1;
 839
 840         l->ave_x=l->sx/l->s1;
 841         l->ave_y=l->sy/l->s1;
 842         l->slope=(l->s1*l->sxy-l->sx*l->sy)/(l->s1*l->sxx-l->sx*l->sx);
 843         l->err_slope=l->yerr/sqrt(sxp);
 844         l->err_ave_y=l->yerr/sqrt(l->s1);
 845
 846         calculate_residuals(l);
 847
 848         return 0;
 849 }
 850
 851 int calculate_residuals(linear_fit_t *l)
 852 {
 853         double x,y;
 854         double res,sum_res2;
 855         int i;
 856
 857         sum_res2=0;
 858         for(i=0;i<l->n;i++){
 859                 x=l->x0+i*l->dx-l->ave_x;
 860                 y=l->y_data[i];
 861
 862                 if(isnan(y))continue;
 863
 864                 res=l->ave_y+l->slope*x-y;
 865
 866                 sum_res2+=res*res;
 867         }
 868         l->S_min=sum_res2/(l->yerr*l->yerr);
 869         l->dof=l->s1-2;
 870
 871         return 0;
 872 }
 873
 874 double linear_fit_func_y(linear_fit_t *l,double x)
 875 {
 876         return l->ave_y+l->slope*(x-l->ave_x);
 877 }
 878
 879
 880 /* printing of scientific numbers (with errors) */
 881
 882 int sci_sprint(char *s,double x,double y)
 883 {
 884         int errsig;
 885         int maxsig;
 886         int sigfigs;
 887         double mantissa;
 888         double error;
 889         double mindigit;
 890
 891         errsig = floor(log10(y));
 892         maxsig = floor(log10(x));
 893         mindigit = pow(10,errsig);
 894
 895         if(maxsig<errsig)maxsig=errsig;
 896
 897         sigfigs = maxsig-errsig+2;
 898
 899         mantissa = x*pow(10,-maxsig);
 900         error = y*pow(10,-errsig+1);
 901
 902         return sprintf(s,"%0.*f(%2.0f)e%d",sigfigs-1,mantissa,error,maxsig);
 903 }
 904
 905 int sci_sprint_alt(char *s,double x,double y)
 906 {
 907         int errsig;
 908         int maxsig;
 909         int sigfigs;
 910         double mantissa;
 911         double error;
 912         double mindigit;
 913
 914         errsig = floor(log10(fabs(y)));
 915         maxsig = floor(log10(fabs(x)));
 916         mindigit = pow(10,errsig);
 917
 918         if(maxsig<errsig)maxsig=errsig;
 919
 920         sigfigs = maxsig-errsig+2;
 921
 922         mantissa = x*pow(10,-maxsig);
 923         error = y*pow(10,-errsig+1);
 924
 925         if(isnan(x)){
 926                 return sprintf(s,"%g",x);
 927         }
 928         if(errsig==1 && maxsig<4 && maxsig>1){
 929                 return sprintf(s,"%0.0f(%2.0f)",x,error);
 930         }
 931         if(maxsig<=0 && maxsig>=-2){
 932                 return sprintf(s,"%0.*f(%2.0f)",sigfigs-1-maxsig,
 933                         mantissa*pow(10,maxsig),error);
 934         }
 935         return sprintf(s,"%0.*f(%2.0f)e%d",sigfigs-1,mantissa,error,maxsig);
 936 }
 937