source: sans/utils/Teabag/sa_rdwr.c @ 838

Last change on this file since 838 was 837, checked in by ajj, 11 years ago

Checkin of Teabag code for TISANE.
File size: 16.7 KB
Line 
1#include <stdio.h>
2#include <math.h>
3#include "sa_data.h"
4#include "sa_futil.h"
5
6/*
7 * I4toI2: Map Integer*4 to Integer*2
8 * Original author: Jim Rhyne
9 * i4toi2 = i4                       if    0 <= i4 <= 32767
10 * i4toi2 = -777                     if    i4 > 2767000
11 * i4toi2 mapped to -13277 to -32768 otherwise
12 *
13 * The mapped values [-776,-1] and [-13276,-778] are not used
14 *
15 */
16short I4toI2 (int i4) {
17  int npw, ipw, ib=10, nd=4, i4max=2767000, i2max = 32767;
18  short i2, error=-777;
19  float r4;
20
21  i2=0;
22  npw=0;
23  if (i4 <= i4max) {
24    if (i4 >= i2max) {
25      r4 = (float) i4;
26      ipw = (int) pow((double)ib,(double)nd);
27      npw = 0;
28      while(r4 >= (ipw - 1.0)) {
29        npw++;
30        r4 /= ib;
31      }
32      i2 = -(short) (r4 + ipw * npw);
33      return i2;
34    }
35    i2 = (short) i4;
36    return i2;
37  }
38  return error;
39}
40
41/*
42 * Adaptation of I2ToI4 routine by Jim Rhyne/Frank Chen
43 * There are better ways of doing this, but for now we should adhere
44 * to the previous method to ensure that it will work.
45 */
46int I2toI4 (short i2) {
47  int i4, npw, ipw, ib=10, nd=4;
48
49  ipw = (int) pow((double) ib,(double) nd);
50  i4 = (int) i2;
51  if (i4 <= -ipw) {
52    npw = -i4 / ipw;
53    i4 = (-i4 % ipw) * ((int) pow((double) ib, (double) npw));
54    return i4;
55  }
56  return i4;
57}
58/*
59 * Read histogram as short int (two byte) and expand to int (four byte),
60 * making adjustments to read with correct byte order.
61 * Little endian : byteorder = 0
62 * Big    endian : byteorder = 1
63 */
64int
65ReadSansHist(FILE *fp, int **intensity, int height, int width, int byteorder) {
66  int x, y;
67  short rawval;
68
69  for(y=height-1;y>=0;y--) {
70    for(x=0;x<width;x++) {
71      if (((ftell(fp) - 514 )%2044) == 0) {
72        if(fgetshort(fp,&rawval,byteorder)) return -1;
73      }
74
75      if(fgetshort(fp,&rawval,byteorder)) return -1;
76      intensity[y][x] = I2toI4(rawval);
77    }
78  }
79  return 0;
80}
81
82int
83WriteSansHist(FILE *fp, int intensity[128][128], int height, int width, int byteorder) 
84{
85  int x, y, offset=0;
86  short rawval, pad=0;
87 
88  for (y=height-1;y>=0;y--) {
89    for(x=0;x<width;x++) {
90      rawval = I4toI2(intensity[y][x]);
91      //printf("%d  %d\n",intensity[y][x],rawval);
92      if ((offset%2044)==0) {
93        if (fputshort(fp,&pad,byteorder)) return -1;
94        offset += 2;
95      }
96      offset += 2;
97      if (fputshort(fp,&rawval,byteorder)) return -1;
98    }
99  }
100  x = ftell(fp);
101  printf("Data written. Offset = %d\n",x);
102  return 0;
103}
104
105int
106ReadSansHdr(FILE *fp, SANSHDR *header, int byteorder) {
107  int i;
108
109  fseek(fp,2,0);
110  if (fgetstring(fp,header->fname,21)) return -1;
111  if (fgetint(fp, &header->run.npre,byteorder)) return -1;
112  if (fgetint(fp, &header->run.ctime,byteorder)) return -1;
113  if (fgetint(fp, &header->run.rtime, byteorder)) return -1;
114  if (fgetint(fp, &header->run.numruns, byteorder)) return -1;
115  if (fgetfloat(fp, &header->run.moncnt, byteorder)) return -1;
116  if (fgetfloat(fp, &header->run.savmon, byteorder)) return -1;
117  if (fgetfloat(fp, &header->run.detcnt, byteorder)) return -1;
118  if (fgetfloat(fp, &header->run.atten, byteorder)) return -1;
119  if (fgetstring(fp,header->run.timdat,20))  return -1;
120  if (fgetstring(fp,header->run.type,3))     return -1;
121  if (fgetstring(fp,header->run.defdir,11))  return -1;
122  fread(&header->run.mode,sizeof(char),1,fp);
123  if (fgetstring(fp,header->run.reserve,8))  return -1;
124  if (fgetstring(fp,header->sample.labl,60)) return -1;
125  if (fgetfloat(fp, &header->sample.trns,byteorder)) return -1;
126  if (fgetfloat(fp, &header->sample.thk,byteorder)) return -1;
127  if (fgetfloat(fp, &header->sample.position,byteorder)) return -1;
128  if (fgetfloat(fp, &header->sample.rotang,byteorder)) return -1;
129  if (fgetint(fp, &header->sample.table,byteorder)) return -1;
130  if (fgetint(fp, &header->sample.holder,byteorder)) return -1;
131  if (fgetint(fp, &header->sample.blank,byteorder)) return -1;
132  if (fgetfloat(fp, &header->sample.temp,byteorder)) return -1;
133  if (fgetfloat(fp, &header->sample.field,byteorder)) return -1;
134  if (fgetint(fp, &header->sample.tctrlr,byteorder)) return -1;
135  if (fgetint(fp, &header->sample.magnet,byteorder)) return -1;
136  if (fgetstring(fp,header->sample.tunits,6)) return -1;
137  if (fgetstring(fp,header->sample.funits,6)) return -1;
138  if (fgetstring(fp,header->det.typ,6)) return -1;
139  for(i=0;i<3;i++) {
140    if (fgetfloat(fp, &header->det.calx[i],byteorder)) return -1;
141  }
142  for(i=0;i<3;i++) {
143    if (fgetfloat(fp, &header->det.caly[i],byteorder)) return -1;
144  }
145  if (fgetint(fp, &header->det.num,byteorder)) return -1;
146  if (fgetint(fp, &header->det.spacer,byteorder)) return -1;
147  if (fgetfloat(fp, &header->det.beamx,byteorder)) return -1;
148  if (fgetfloat(fp, &header->det.beamy,byteorder)) return -1;
149  if (fgetfloat(fp, &header->det.dis,byteorder)) return -1;
150  if (fgetfloat(fp, &header->det.ang,byteorder)) return -1;
151  if (fgetfloat(fp, &header->det.siz,byteorder)) return -1;
152  if (fgetfloat(fp, &header->det.bstop,byteorder)) return -1;
153  if (fgetfloat(fp, &header->det.blank,byteorder)) return -1;
154  if (fgetfloat(fp, &header->resolution.ap1,byteorder)) return -1;
155  if (fgetfloat(fp, &header->resolution.ap2,byteorder)) return -1;
156  if (fgetfloat(fp, &header->resolution.ap12dis,byteorder)) return -1;
157  if (fgetfloat(fp, &header->resolution.lmda,byteorder)) return -1;
158  if (fgetfloat(fp, &header->resolution.dlmda,byteorder)) return -1;
159  if (fgetfloat(fp, &header->resolution.save,byteorder)) return -1;
160  if (fgetbool(fp, &header->tslice.slicing,byteorder)) return -1;
161  if (fgetint(fp, &header->tslice.multfact,byteorder)) return -1;
162  if (fgetint(fp, &header->tslice.ltslice,byteorder)) return -1;
163  if (fgetbool(fp, &header->temp.printemp,byteorder)) return -1;
164  if (fgetfloat(fp, &header->temp.hold,byteorder)) return -1;
165  if (fgetfloat(fp, &header->temp.err,byteorder)) return -1;
166  if (fgetfloat(fp, &header->temp.blank,byteorder)) return -1;
167  if (fgetint(fp, &header->temp.extra,byteorder)) return -1;
168  if (fgetint(fp, &header->temp.reserve,byteorder)) return -1;
169  if (fgetbool(fp, &header->magnet.printmag,byteorder)) return -1;
170  if (fgetbool(fp, &header->magnet.sensor,byteorder)) return -1;
171  if (fgetfloat(fp, &header->magnet.current,byteorder)) return -1;
172  if (fgetfloat(fp, &header->magnet.conv,byteorder)) return -1;
173  if (fgetfloat(fp, &header->magnet.fieldlast,byteorder)) return -1;
174  if (fgetfloat(fp, &header->magnet.blank,byteorder)) return -1;
175  if (fgetfloat(fp, &header->magnet.spacer,byteorder)) return -1;
176  if (fgetfloat(fp, &header->bmstp.xpos,byteorder)) return -1;
177  if (fgetfloat(fp, &header->bmstp.ypos,byteorder)) return -1;
178  if (fgetint(fp, &header->params.blank1,byteorder)) return -1;
179  if (fgetint(fp, &header->params.blank2,byteorder)) return -1;
180  if (fgetint(fp, &header->params.blank3,byteorder)) return -1;
181  if (fgetfloat(fp, &header->params.trnscnt,byteorder)) return -1;
182  if (fgetfloat(fp, &header->params.extra1,byteorder)) return -1;
183  if (fgetfloat(fp, &header->params.extra2,byteorder)) return -1;
184  if (fgetfloat(fp, &header->params.extra3,byteorder)) return -1;
185  if (fgetstring(fp,header->params.reserve,42)) return -1;
186  if (fgetbool(fp, &header->voltage.printvolt,byteorder)) return -1;
187  if (fgetfloat(fp, &header->voltage.volts,byteorder)) return -1;
188  if (fgetfloat(fp, &header->voltage.blank,byteorder)) return -1;
189  if (fgetint(fp, &header->voltage.spacer,byteorder)) return -1;
190  if (fgetbool(fp, &header->polarization.printpol,byteorder)) return -1;
191  if (fgetbool(fp, &header->polarization.flipper,byteorder)) return -1;
192  if (fgetfloat(fp, &header->polarization.horiz,byteorder)) return -1;
193  if (fgetfloat(fp, &header->polarization.vert,byteorder)) return -1;
194  for(i=0;i<2;i++) {
195    if (fgetint(fp, &header->analysis.rows[i],byteorder)) return -1;
196  }
197  for(i=0;i<2;i++) {
198    if (fgetint(fp, &header->analysis.cols[i],byteorder)) return -1;
199  }
200  if (fgetfloat(fp, &header->analysis.factor,byteorder)) return -1;
201  if (fgetfloat(fp, &header->analysis.qmin,byteorder)) return -1;
202  if (fgetfloat(fp, &header->analysis.qmax,byteorder)) return -1;
203  if (fgetfloat(fp, &header->analysis.imin,byteorder)) return -1;
204  if (fgetfloat(fp, &header->analysis.imax,byteorder)) return -1;
205  i = ftell(fp);
206  printf("Header read: offset = %d\n",i);
207  return 0;
208}
209
210int
211WriteSansHdr(FILE *fp, SANSHDR *header, int byteorder) {
212  int i;
213  //short spad=0;
214
215  fwrite(fp,sizeof(short),1,fp);
216  if (fputstring(fp,header->fname,21)) return -1;
217  if (fputint(fp, &header->run.npre,byteorder)) return -1;
218  if (fputint(fp, &header->run.ctime,byteorder)) return -1;
219  if (fputint(fp, &header->run.rtime, byteorder)) return -1;
220  if (fputint(fp, &header->run.numruns, byteorder)) return -1;
221  if (fputfloat(fp, &header->run.moncnt, byteorder)) return -1;
222  if (fputfloat(fp, &header->run.savmon, byteorder)) return -1;
223  if (fputfloat(fp, &header->run.detcnt, byteorder)) return -1;
224  if (fputfloat(fp, &header->run.atten, byteorder)) return -1;
225  if (fputstring(fp,header->run.timdat,20))  return -1;
226  if (fputstring(fp,header->run.type,3))     return -1;
227  if (fputstring(fp,header->run.defdir,11))  return -1;
228  fwrite(&header->run.mode,sizeof(char),1,fp);
229  if (fputstring(fp,header->run.reserve,8))  return -1;
230  if (fputstring(fp,header->sample.labl,60)) return -1;
231  if (fputfloat(fp, &header->sample.trns,byteorder)) return -1;
232  if (fputfloat(fp, &header->sample.thk,byteorder)) return -1;
233  if (fputfloat(fp, &header->sample.position,byteorder)) return -1;
234  if (fputfloat(fp, &header->sample.rotang,byteorder)) return -1;
235  if (fputint(fp, &header->sample.table,byteorder)) return -1;
236  if (fputint(fp, &header->sample.holder,byteorder)) return -1;
237  if (fputint(fp, &header->sample.blank,byteorder)) return -1;
238  if (fputfloat(fp, &header->sample.temp,byteorder)) return -1;
239  if (fputfloat(fp, &header->sample.field,byteorder)) return -1;
240  if (fputint(fp, &header->sample.tctrlr,byteorder)) return -1;
241  if (fputint(fp, &header->sample.magnet,byteorder)) return -1;
242  if (fputstring(fp,header->sample.tunits,6)) return -1;
243  if (fputstring(fp,header->sample.funits,6)) return -1;
244  if (fputstring(fp,header->det.typ,6)) return -1;
245  for(i=0;i<3;i++) {
246    if (fputfloat(fp, &header->det.calx[i],byteorder)) return -1;
247  }
248  for(i=0;i<3;i++) {
249    if (fputfloat(fp, &header->det.caly[i],byteorder)) return -1;
250  }
251  if (fputint(fp, &header->det.num,byteorder)) return -1;
252  if (fputint(fp, &header->det.spacer,byteorder)) return -1;
253  if (fputfloat(fp, &header->det.beamx,byteorder)) return -1;
254  if (fputfloat(fp, &header->det.beamy,byteorder)) return -1;
255  if (fputfloat(fp, &header->det.dis,byteorder)) return -1;
256  if (fputfloat(fp, &header->det.ang,byteorder)) return -1;
257  if (fputfloat(fp, &header->det.siz,byteorder)) return -1;
258  if (fputfloat(fp, &header->det.bstop,byteorder)) return -1;
259  if (fputfloat(fp, &header->det.blank,byteorder)) return -1;
260  if (fputfloat(fp, &header->resolution.ap1,byteorder)) return -1;
261  if (fputfloat(fp, &header->resolution.ap2,byteorder)) return -1;
262  if (fputfloat(fp, &header->resolution.ap12dis,byteorder)) return -1;
263  if (fputfloat(fp, &header->resolution.lmda,byteorder)) return -1;
264  if (fputfloat(fp, &header->resolution.dlmda,byteorder)) return -1;
265  if (fputfloat(fp, &header->resolution.save,byteorder)) return -1;
266  if (fputbool(fp, &header->tslice.slicing,byteorder)) return -1;
267  if (fputint(fp, &header->tslice.multfact,byteorder)) return -1;
268  if (fputint(fp, &header->tslice.ltslice,byteorder)) return -1;
269  if (fputbool(fp, &header->temp.printemp,byteorder)) return -1;
270  if (fputfloat(fp, &header->temp.hold,byteorder)) return -1;
271  if (fputfloat(fp, &header->temp.err,byteorder)) return -1;
272  if (fputfloat(fp, &header->temp.blank,byteorder)) return -1;
273  if (fputint(fp, &header->temp.extra,byteorder)) return -1;
274  if (fputint(fp, &header->temp.reserve,byteorder)) return -1;
275  if (fputbool(fp, &header->magnet.printmag,byteorder)) return -1;
276  if (fputbool(fp, &header->magnet.sensor,byteorder)) return -1;
277  if (fputfloat(fp, &header->magnet.current,byteorder)) return -1;
278  if (fputfloat(fp, &header->magnet.conv,byteorder)) return -1;
279  if (fputfloat(fp, &header->magnet.fieldlast,byteorder)) return -1;
280  if (fputfloat(fp, &header->magnet.blank,byteorder)) return -1;
281  if (fputfloat(fp, &header->magnet.spacer,byteorder)) return -1;
282  if (fputfloat(fp, &header->bmstp.xpos,byteorder)) return -1;
283  if (fputfloat(fp, &header->bmstp.ypos,byteorder)) return -1;
284  if (fputint(fp, &header->params.blank1,byteorder)) return -1;
285  if (fputint(fp, &header->params.blank2,byteorder)) return -1;
286  if (fputint(fp, &header->params.blank3,byteorder)) return -1;
287  if (fputfloat(fp, &header->params.trnscnt,byteorder)) return -1;
288  if (fputfloat(fp, &header->params.extra1,byteorder)) return -1;
289  if (fputfloat(fp, &header->params.extra2,byteorder)) return -1;
290  if (fputfloat(fp, &header->params.extra3,byteorder)) return -1;
291  if (fputstring(fp,header->params.reserve,42)) return -1;
292  if (fputbool(fp, &header->voltage.printvolt,byteorder)) return -1;
293  if (fputfloat(fp, &header->voltage.volts,byteorder)) return -1;
294  if (fputfloat(fp, &header->voltage.blank,byteorder)) return -1;
295  if (fputint(fp, &header->voltage.spacer,byteorder)) return -1;
296  if (fputbool(fp, &header->polarization.printpol,byteorder)) return -1;
297  if (fputbool(fp, &header->polarization.flipper,byteorder)) return -1;
298  if (fputfloat(fp, &header->polarization.horiz,byteorder)) return -1;
299  if (fputfloat(fp, &header->polarization.vert,byteorder)) return -1;
300  for(i=0;i<2;i++) {
301    if (fputint(fp, &header->analysis.rows[i],byteorder)) return -1;
302  }
303  for(i=0;i<2;i++) {
304    if (fputint(fp, &header->analysis.cols[i],byteorder)) return -1;
305  }
306  if (fputfloat(fp, &header->analysis.factor,byteorder)) return -1;
307  if (fputfloat(fp, &header->analysis.qmin,byteorder)) return -1;
308  if (fputfloat(fp, &header->analysis.qmax,byteorder)) return -1;
309  if (fputfloat(fp, &header->analysis.imin,byteorder)) return -1;
310  if (fputfloat(fp, &header->analysis.imax,byteorder)) return -1;
311  i = ftell(fp);
312  printf("Header written: offset = %d\n",i);
313  return 0;
314}
315
316int
317DumpSansHdr(SANSHDR *header) {
318  int i;
319
320  printf("RUN: \n");
321  printf("  npre   = %d\n", header->run.npre);
322  printf("  ctime  = %d\n", header->run.ctime);
323  printf("  rtime  = %d\n", header->run.rtime);
324  printf("  numruns= %d\n", header->run.numruns);
325  printf("  moncnt = %f\n", header->run.moncnt);
326  printf("  savmon = %f\n", header->run.savmon);
327  printf("  detcnt = %f\n", header->run.detcnt);
328  printf("  atten  = %f\n", header->run.atten);
329  printf("  timdat = '%s'\n",header->run.timdat);
330  printf("  type   = '%s'\n",header->run.type);
331  printf("  defdir = '%s'\n",header->run.defdir);
332
333  printf("DET: \n");
334  printf("  typ    = '%s'\n",header->det.typ);
335  for (i=0;i<3;i++) printf("  calx[%d] = %f\n",i+1,header->det.calx[i]);
336  for (i=0;i<3;i++) printf("  caly[%d] = %f\n",i+1,header->det.caly[i]);
337  printf("  beamx  = %f\n",header->det.beamx);
338  printf("  beamy  = %f\n",header->det.beamy);
339  printf("  dis    = %f\n",header->det.dis);
340  printf("  ang    = %f\n",header->det.ang);
341  printf("  siz    = %f\n",header->det.siz);
342  printf("  bstop  = %f\n",header->det.bstop);
343
344  printf("RESOLUTION: \n");
345  printf("  ap1    = %f\n",header->resolution.ap1);
346  printf("  ap2    = %f\n",header->resolution.ap2);
347  printf("  ap12dis= %f\n",header->resolution.ap12dis);
348  printf("  lmda   = %f\n",header->resolution.lmda);
349  printf("  dlmda  = %f\n",header->resolution.dlmda);
350  printf("  save   = %f\n",header->resolution.save);
351
352  printf("BEAMSTOP: \n");
353  printf("  xpos   = %f\n",header->bmstp.xpos);
354  printf("  ypos   = %f\n",header->bmstp.ypos);
355
356  return 0;
357}
358
359int InitSansHdr(SANSHDR * header) 
360{
361  memset(header, 0,sizeof(SANSHDR));
362  strcpy(header->fname,"default.sans");
363  header->run.npre = 0;
364  header->run.ctime = 0;
365  return 0;
366}
367
368#ifndef PI
369#define PI 3.141529
370#endif
371
372double gaussian(double x, double x0, double amp, double fwhm)
373{
374  double xprime; 
375
376  xprime = 2*(x - x0)/fwhm;
377  return amp * exp( -1 * (xprime)*(xprime));
378}
379
380double lorentzian(double x, double x0, double amp, double fwhm)
381{
382  double xprime;
383  xprime = 2*(x - x0)/fwhm;
384  return amp / (1 + xprime * xprime);
385}
386
387/* This could stand some improvement */
388int FakeSansHist(int hist[128][128], double distance, double dim, double lambda)
389{
390  double peak = 1000;
391  double ki,q,qmin;
392  int x0=63, y0=63, i, j;
393  //double x, y;
394  double len, r, theta, intensity;
395 
396  ki = 2*PI/lambda;
397  len = dim/128; /* Length of one interval */
398
399  printf("lambda = %f\n",lambda);
400  printf("ki = %f\n",ki);
401  r = len * sqrt(2) * 64;
402  theta = atan2(r,distance)/2;
403  q = 2 * ki * sin(theta);
404  printf("maxq = %f\n",q);
405
406  theta = atan2(len,distance)/2;
407  qmin = 2 * ki * sin(theta);
408  printf("minq = %f\n",qmin);
409
410  for(i=0;i<128;i++) {
411    for (j=0;j<128;j++) {
412      r = len * sqrt((i-x0)*(i-x0) + (j-y0)*(j-y0));
413      theta = atan2(r,distance)/2;
414      q = 2 * ki * sin(theta);
415      if (q == 0) { intensity = peak; } else {intensity = peak* (qmin*qmin)/(q*q);}
416      hist[i][j] = (int) intensity;
417    }
418  }
419  return 0;
420}
Note: See TracBrowser for help on using the repository browser.