Ignore:
Timestamp:
Oct 26, 2010 3:32:14 PM (12 years ago)
Author:
srkline
Message:

Changes to DebyeSpheres? to convert everything to double, rather than float.

Changes to MonteCarlo? to use a local definition of round(), since it's not a standard function in Visual Studio's math.h. Easier to just write my own and give it an odd name - MC_round()

Location:
sans/XOP_Dev/MonteCarlo
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • sans/XOP_Dev/MonteCarlo/DebyeSpheres.c

    r623 r758  
    2222    using Debye's method.  
    2323 
    24         Warning: 
    25                 The call to WaveData() below returns a pointer to the middle 
    26                 of an unlocked Macintosh handle. In the unlikely event that your 
    27                 calculations could cause memory to move, you should copy the coefficient 
    28                 values to local variables or an array before such operations. 
     24 everything was previously declared as float, some inputs were as double. 
     25  
     26 -- to avoid confusion, calculate everything in double, pass everything in as double 
     27  -- this is more precision than is necessary, but avoids strange results from type casting 
     28     behind the scenes 
    2929*/ 
    3030int 
    3131DebyeSpheresX(AltiParamsPtr p) 
    3232{ 
    33         float qv;                               // input q-value 
    34         float ival;                             //output intensity value 
    35         float *xv,*yv,*zv,*rv;          //pointers to input xyz-rho coordinates 
     33        double qv;                              // input q-value 
     34        double ival;                            //output intensity value 
     35        double *xv,*yv,*zv,*rv;         //pointers to input xyz-rho coordinates 
    3636        int i,j; 
    3737    int npt; 
    38         float rval,grid,vol,fQR,dum,dij; 
     38        double rval,grid,vol,fQR,dum,dij; 
    3939 
    4040         
    … …  
    6060 
    6161    //check to see that all are float, not double 
     62        /* 
    6263        if(WaveType(p->rhowavH) != NT_FP32 ) { 
    6364        SetNaN64(&p->result); 
    … …  
    7677        return kExpectedNT_FP32; 
    7778    } 
    78  
     79        */ 
    7980         
    80         // convert the input to float. Do all calculations as float. 
    81  
    82     rval = (float) p->Rprimary;         // primary sphere radius 
    83     grid = (float) p->grid;                     // calling program should set this to 0.62*Rprimary 
    84         qv = (float) p->qval; 
     81    //check to see that all are double 
     82        if(WaveType(p->rhowavH) != NT_FP64 ) { 
     83        SetNaN64(&p->result); 
     84        return kExpectedNT_FP64; 
     85    } 
     86    if(WaveType(p->zwavH) != NT_FP64 ) { 
     87        SetNaN64(&p->result); 
     88        return kExpectedNT_FP64; 
     89    } 
     90    if(WaveType(p->ywavH) != NT_FP64 ) { 
     91        SetNaN64(&p->result); 
     92        return kExpectedNT_FP64; 
     93    } 
     94    if(WaveType(p->xwavH) != NT_FP64 ) { 
     95        SetNaN64(&p->result); 
     96        return kExpectedNT_FP64; 
     97    } 
     98         
     99         
     100        // (NO) -- convert the input to float. Do all calculations as float. 
     101        // do everything in double. no reason to do in float, except for the previous Altivec incarnation 
     102         
     103    rval = p->Rprimary;         // primary sphere radius 
     104    grid = p->grid;                     // calling program should set this to 0.62*Rprimary 
     105        qv = p->qval; 
    85106 
    86107// 
    … …  
    109130        } 
    110131 
    111     p->result= (DOUBLE) ival; 
     132    p->result = ival; 
    112133         
    113134    return 0; 
    … …  
    116137 
    117138 
    118 float PhiQR(float qval, float rval) 
     139double PhiQR(double qval, double rval) 
    119140{ 
    120         float retval,qr; 
     141        double retval,qr; 
    121142         
    122143        qr = qval*rval; 
    123         retval = (float) ( 3*(sin(qr) - qr*cos(qr))/qr/qr/qr ); 
     144        retval = ( 3*(sin(qr) - qr*cos(qr))/qr/qr/qr ); 
    124145        return(retval); 
    125146} 
    126147 
    127 float XYZDistance(float x1, float x2,float y1, float y2,float z1, float z2) 
     148double XYZDistance(double x1, double x2,double y1, double y2,double z1, double z2) 
    128149{ 
    129         float retval,dx,dy,dz; 
     150        double retval,dx,dy,dz; 
    130151         
    131152        dx = (x1-x2); 
    132153        dy = (y1-y2); 
    133154        dz = (z1-z2); 
    134         retval = (float) sqrt( dx*dx + dy*dy + dz*dz ); 
     155        retval = sqrt( dx*dx + dy*dy + dz*dz ); 
    135156        return(retval); 
    136157} 

  • sans/XOP_Dev/MonteCarlo/DebyeSpheres.h

    r674 r758  
    1212// result is the last parameter, always. 
    1313typedef struct AltiParams { 
    14     DOUBLE       grid;          // effective c-to-c distance between spheres = 0.62*Rprimary 
    15     DOUBLE   Rprimary;  //primary sphere radius 
     14    double       grid;          // effective c-to-c distance between spheres = 0.62*Rprimary 
     15    double   Rprimary;  //primary sphere radius 
    1616        waveHndl rhowavH;       // rho at xyz!!! 
    1717        waveHndl zwavH; // z coordinate. ALL are expected to be SP waves 
    1818        waveHndl ywavH; // y coordinate. 
    1919        waveHndl xwavH; // x coordinate. 
    20     DOUBLE qval;        // q-value. 
     20    double qval;        // q-value. 
    2121        void* tp;                       //unused void for threadsafe functions 
    22         DOUBLE result; 
     22        double result; 
    2323}AltiParams, *AltiParamsPtr;    
    2424  
    2525  
    26 float PhiQR(float qval, float rval); 
    27 float XYZDistance(float x1, float x2,float y1, float y2,float z1, float z2); 
     26double PhiQR(double qval, double rval); 
     27double XYZDistance(double x1, double x2,double y1, double y2,double z1, double z2); 
    2828int DebyeSpheresX(AltiParamsPtr p); 

  • sans/XOP_Dev/MonteCarlo/MonteCarlo.h

    r623 r758  
    5454int Monte_SANSX3(MC_ParamsPtr p); 
    5555int Monte_SANSX4(MC_ParamsPtr p); 
     56 
     57long MC_round(double x); 
    5658int FindPixel(double testQ, double testPhi, double lam, double sdd, double pixSize, double xCtr, double yCtr, long *xPixel, long *yPixel); 
    5759int NewDirection(double *vx, double *vy, double *vz, double theta, double phi); 

  • sans/XOP_Dev/MonteCarlo/MonteCarlo.r

    r673 r758  
    6868                { 
    6969                                                NT_FP64,                        /* single precision wave (q-wave) */ 
    70                         NT_FP32 + WAVE_TYPE,                            /* single precision wave (x coordinates) */ 
    71                         NT_FP32 + WAVE_TYPE,                            /* single precision wave (y coordinates) */ 
    72                         NT_FP32 + WAVE_TYPE,                            /* single precision wave (z coordinates) */ 
    73                         NT_FP32 + WAVE_TYPE,                            /* single precision wave (rho at xyz) */ 
     70                        NT_FP64 + WAVE_TYPE,                            /* single precision wave (x coordinates) */ 
     71                        NT_FP64 + WAVE_TYPE,                            /* single precision wave (y coordinates) */ 
     72                        NT_FP64 + WAVE_TYPE,                            /* single precision wave (z coordinates) */ 
     73                        NT_FP64 + WAVE_TYPE,                            /* single precision wave (rho at xyz) */ 
    7474                        NT_FP64,                                // Rprimary, the primary sphere radius 
    7575                        NT_FP64,                                // grid, should be passed as 0.62*Rprimary 

  • sans/XOP_Dev/MonteCarlo/MonteCarlo_Main.c

    r711 r758  
    1616#include "DebyeSpheres.h" 
    1717 
     18 
     19// this function is here simply because MS visual studio does not contain lround() or round() in math.h 
     20// rounds away from zero 
     21// -- only used in FindPixel 
     22long MC_round(double x) 
     23{ 
     24        if(x == 0)      { 
     25                return (long)(0); 
     26        } 
     27        if(x > 0)       { 
     28                return (long)(x + 0.5f); 
     29        } else {                // x < 0 
     30                return (long)(x - 0.5f); 
     31        } 
     32} 
     33 
     34 
    1835int 
    1936FindPixel(double testQ, double testPhi, double lam, double sdd, 
    … …  
    2946        theta = 2.0*asin(qy*lam/4.0/pi); 
    3047        dy = sdd*tan(theta); 
    31         *yPixel = lround(yCtr + dy/pixSize);            //corrected 7/2010 to round away from zero, to avoid 2x counts in row 0 and column 0 
     48//      *yPixel = lround(yCtr + dy/pixSize);            //corrected 7/2010 to round away from zero, to avoid 2x counts in row 0 and column 0 
     49        *yPixel = MC_round(yCtr + dy/pixSize);          //corrected 7/2010 to round away from zero, to avoid 2x counts in row 0 and column 0 
    3250         
    3351        theta = 2.0*asin(qx*lam/4.0/pi); 
    3452        dx = sdd*tan(theta); 
    35         *xPixel = lround(xCtr + dx/pixSize); 
     53//      *xPixel = lround(xCtr + dx/pixSize); 
     54        *xPixel = MC_round(xCtr + dx/pixSize); 
    3655         
    3756        //if on detector, return xPix and yPix values, otherwise -1 
Note: See TracChangeset for help on using the changeset viewer.