Changeset 6e93a02 in sasview for sansmodels/src/libigor/libSphere.c


Ignore:
Timestamp:
Mar 30, 2010 7:53:47 PM (14 years ago)
Author:
Jae Cho <jhjcho@…>
Branches:
master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
Children:
f10063e
Parents:
64017a8
Message:

updated libigor files from NIST svn

File:
1 edited

Legend:

Unmodified
Added
Removed

  • sansmodels/src/libigor/libSphere.c

    r975ec8e r6e93a02  
    1313SphereForm(double dp[], double q) 
    1414{ 
    15         double scale,radius,delrho,bkg;         //my local names 
     15        double scale,radius,delrho,bkg,sldSph,sldSolv;          //my local names 
    1616        double bes,f,vol,f2,pi,qr; 
    1717 
    … …  
    1919        scale = dp[0]; 
    2020        radius = dp[1]; 
    21         delrho = dp[2]; 
    22         bkg = dp[3]; 
    23  
     21        sldSph = dp[2]; 
     22        sldSolv = dp[3]; 
     23        bkg = dp[4]; 
     24         
     25        delrho = sldSph - sldSolv; 
     26        //handle qr==0 separately 
    2427        qr = q*radius; 
    25         if(qr == 0){ 
     28        if(qr == 0.0){ 
    2629                bes = 1.0; 
    2730        }else{ 
    2831                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
    2932        } 
    30  
     33         
    3134        vol = 4.0*pi/3.0*radius*radius*radius; 
    3235        f = vol*bes*delrho;             // [=] A-1 
    3336                                                        // normalize to single particle volume, convert to 1/cm 
    3437        f2 = f * f / vol * 1.0e8;               // [=] 1/cm 
    35  
     38         
    3639        return(scale*f2+bkg);   //scale, and add in the background 
    3740} 
    … …  
    4447        double scale,rcore,thick,rhocore,rhoshel,rhosolv,bkg;           //my local names 
    4548        double bes,f,vol,qr,contr,f2; 
    46  
     49         
    4750        pi = 4.0*atan(1.0); 
    4851        x=q; 
    49  
     52         
    5053        scale = dp[0]; 
    5154        rcore = dp[1]; 
    … …  
    5861        qr=x*rcore; 
    5962        contr = rhocore-rhoshel; 
    60         if(qr == 0){ 
     63        if(qr == 0.0){ 
    6164                bes = 1.0; 
    6265        }else{ 
    … …  
    6871        qr=x*(rcore+thick); 
    6972        contr = rhoshel-rhosolv; 
    70         if(qr == 0){ 
     73        if(qr == 0.0){ 
    7174                bes = 1.0; 
    7275        }else{ 
    … …  
    7881        // normalize to particle volume and rescale from [A-1] to [cm-1] 
    7982        f2 = f*f/vol*1.0e8; 
    80  
     83         
    8184        //scale if desired 
    8285        f2 *= scale; 
    8386        // then add in the background 
    8487        f2 += bkg; 
    85  
     88         
    8689        return(f2); 
    8790} 
    … …  
    97100        pi = 4.0*atan(1.0); 
    98101        x= q; 
    99  
     102         
    100103        scale = dp[0]; 
    101104        rcore = dp[1]; 
    … …  
    118121        qr=x*(rcore+thick); 
    119122        contr = rhoshel-rhosolv; 
    120         if(qr == 0){ 
     123        if(qr == 0.0){ 
    121124                bes = 1.0; 
    122125        }else{ 
    … …  
    130133        vol = 4.0*pi/3.0*(pow((rcore+thick),3)-pow(rcore,3)); 
    131134        f2 = f*f/vol*1.0e8; 
    132  
     135         
    133136        //scale if desired 
    134137        f2 *= scale; 
    135138        // then add in the background 
    136139        f2 += bkg; 
    137  
     140         
    138141        return(f2); 
    139142} 
    … …  
    152155        double zp3,vpoly; 
    153156        double s2y, arg1, arg2, arg3, drh1, drh2; 
    154  
     157         
    155158        pi = 4.0*atan(1.0); 
    156159        qq= q; 
    … …  
    162165        drho2 = dp[5]-dp[6];            //shell-solvent 
    163166        bkg = dp[7]; 
    164  
    165         zz = (1.0/sig)*(1.0/sig) - 1.0; 
    166  
     167         
     168        zz = (1.0/sig)*(1.0/sig) - 1.0;    
     169         
    167170        h=qq; 
    168  
     171         
    169172        drh1 = drho1; 
    170173        drh2 = drho2; 
    … …  
    174177        zp3 = zz + 3.; 
    175178        vpoly = 4*pi/3*zp3*zp2/zp1/zp1*pow((corrad+del),3); 
    176  
    177  
     179         
     180         
    178181        // remember that h is the passed in value of q for the calculation 
    179182        y = h *del; 
    … …  
    196199        c8 = c4 - .5 + g * cy - g * .5 * g * (y * y + 1.) * c2y; 
    197200        c9 = g * sy * (1. - g * cy); 
    198  
     201         
    199202        tb = log(zp1 * zp1 / (zp1 * zp1 + x * 4. * x)); 
    200203        tb1 = exp(zp1 * .5 * tb); 
    201204        tb2 = exp(zp2 * .5 * tb); 
    202205        tb3 = exp(zp3 * .5 * tb); 
    203  
     206         
    204207        t1 = c1 + c2 * x + c3 * x * x * zp2 / zp1; 
    205208        t2 = tb1 * (c4 * cos(arg1) + c7 * sin(arg1)); 
    … …  
    216219        // then add in the background 
    217220        form += bkg; 
    218  
     221         
    219222        return(form); 
    220223} 
    … …  
    234237        double capl,capl1,capmu,capmu1,r3,pq; 
    235238        double ka2,r1,heff; 
    236         double hh,k; 
    237  
     239        double hh,k,slds,sld; 
     240         
    238241        pi = 4.0*atan(1.0); 
    239242        k= q; 
    240  
     243         
    241244        rad = dp[0];            // radius (A) 
    242245        z2 = dp[1];             //polydispersity (0<z2<1) 
    243246        phi = dp[2];            // volume fraction (0<phi<1) 
    244         cont = dp[3]*1.0e4;             // contrast (odd units) 
    245         bkg = dp[4]; 
     247        slds = dp[3]; 
     248        sld = dp[4]; 
     249        cont = (slds - sld)*1.0e4;              // contrast (odd units) 
     250        bkg = dp[5]; 
    246251        sigma = 2*rad; 
    247  
     252         
    248253        zz=1.0/(z2*z2)-1.0; 
    249254        bb = sigma/(zz+1.0); 
    250255        cc = zz+1.0; 
    251  
     256         
    252257        //*c   Compute the number density by <r-cubed>, not <r> cubed*/ 
    253258        r1 = sigma/2.0; 
    … …  
    281286        d5=pow((pi/capd),2)*((rho/k)*(chi-1.0)+0.5*k*t2); 
    282287        d6=pow((pi/capd),2)*(rho/k)*(chi2-s2); 
    283  
     288         
    284289        e1=pow((pi/capd),2)*pow((rho/k/k),2)*((chi-1.0)*(chi2-s2)-(chi1-s1)*(chi1-s1)-(k*s1-pp)*(k*s3-p2)+pow((k*s2-p1),2)); 
    285290        ee=1.0-(2.0*pi/capd)*(1.0+0.5*pi*t3/capd)*(rho/k/k/k)*(k*s1-pp)-(2.0*pi/capd)*rho/k/k*((chi1-s1)+(0.25*pi*t2/capd)*(chi2-s2))-e1; 
    286291        y1=pow((pi/capd),2)*pow((rho/k/k),2)*((k*s1-pp)*(chi2-s2)-2.0*(k*s2-p1)*(chi1-s1)+(k*s3-p2)*(chi-1.0)); 
    287         yy = (2.0*pi/capd)*(1.0+0.5*pi*t3/capd)*(rho/k/k/k)*(chi+0.5*k*k*s2-1.0)-(2.0*pi*rho/capd/k/k)*(k*s2-p1+(0.25*pi*t2/capd)*(k*s3-p2))-y1; 
    288  
     292        yy = (2.0*pi/capd)*(1.0+0.5*pi*t3/capd)*(rho/k/k/k)*(chi+0.5*k*k*s2-1.0)-(2.0*pi*rho/capd/k/k)*(k*s2-p1+(0.25*pi*t2/capd)*(k*s3-p2))-y1;     
     293         
    289294        capl=2.0*pi*cont*rho/k/k/k*(pp-0.5*k*(s1+chi1)); 
    290295        capl1=2.0*pi*cont*rho/k/k/k*(p1-0.5*k*(s2+chi2)); 
    291296        capmu=2.0*pi*cont*rho/k/k/k*(1.0-chi-0.5*k*p1); 
    292297        capmu1=2.0*pi*cont*rho/k/k/k*(s1-chi1-0.5*k*p2); 
    293  
     298         
    294299        h1=capl*(capl*(yy*d1-ee*d6)+capl1*(yy*d2-ee*d4)+capmu*(ee*d1+yy*d6)+capmu1*(ee*d2+yy*d4)); 
    295300        h2=capl1*(capl*(yy*d2-ee*d4)+capl1*(yy*d3-ee*d5)+capmu*(ee*d2+yy*d4)+capmu1*(ee*d3+yy*d5)); 
    296301        h3=capmu*(capl*(ee*d1+yy*d6)+capl1*(ee*d2+yy*d4)+capmu*(ee*d6-yy*d1)+capmu1*(ee*d4-yy*d2)); 
    297302        h4=capmu1*(capl*(ee*d2+yy*d4)+capl1*(ee*d3+yy*d5)+capmu*(ee*d4-yy*d2)+capmu1*(ee*d5-yy*d3)); 
    298  
     303         
    299304        //*  This part computes the second integral in equation (1) of the paper.*/ 
    300  
     305         
    301306        hint1 = -2.0*(h1+h2+h3+h4)/(k*k*k*(ee*ee+yy*yy)); 
    302  
     307         
    303308        //*  This part computes the first integral in equation (1).  It also 
    304309        // generates the KC approximated effective structure factor.*/ 
    305  
     310         
    306311        pq=4.0*pi*cont*(sin(k*sigma/2.0)-0.5*k*sigma*cos(k*sigma/2.0)); 
    307312        hint2=8.0*pi*pi*rho*cont*cont/(k*k*k*k*k*k)*(1.0-chi-k*p1+0.25*k*k*(s2+chi2)); 
    308  
     313         
    309314        ka=k*(sigma/2.0); 
    310315        // 
    … …  
    314319        ka2=ka*ka; 
    315320        //* 
    316         //  heff is PY analytical solution for intensity divided by the 
     321        //  heff is PY analytical solution for intensity divided by the  
    317322        //   form factor.  happ is the KC approximated effective S(q) 
    318  
     323          
    319324         //******************* 
    320325         //  add in the background then return the intensity value 
    321  
     326          
    322327         return(hh+bkg);        //scale, and add in the background 
    323328} 
    … …  
    332337        double scale,ravg,pd,bkg,rho,rhos;              //my local names 
    333338        double scale2,ravg2,pd2,rho2;           //my local names 
    334  
     339         
    335340        x= q; 
    336  
     341         
    337342        scale = dp[0]; 
    338343        ravg = dp[1]; 
    … …  
    345350        rhos = dp[8]; 
    346351        bkg = dp[9]; 
    347  
     352         
    348353        pq = SchulzSphere_Fn( scale,  ravg,  pd,  rho,  rhos, x); 
    349354        pq += SchulzSphere_Fn( scale2,  ravg2,  pd2,  rho2,  rhos, x); 
    350355        // add in the background 
    351356        pq += bkg; 
    352  
     357         
    353358        return (pq); 
    354359} 
    … …  
    361366        double x,pq; 
    362367        double scale,ravg,pd,bkg,rho,rhos;              //my local names 
    363  
     368         
    364369        x= q; 
    365  
     370         
    366371        scale = dp[0]; 
    367372        ravg = dp[1]; 
    … …  
    373378        // add in the background 
    374379        pq += bkg; 
    375  
     380         
    376381        return(pq); 
    377382} 
    … …  
    384389        double aa,at1,at2,rt1,rt2,rt3,t1,t2,t3; 
    385390        double v1,v2,v3,g1,pq,pi,delrho,zz; 
    386  
     391        double i_zero,Rg2,zp8; 
     392         
    387393        pi = 4.0*atan(1.0); 
    388394        delrho = rho-rhos; 
    389         zz = (1.0/pd)*(1.0/pd) - 1.0; 
    390  
     395        zz = (1.0/pd)*(1.0/pd) - 1.0;    
     396         
    391397        zp1 = zz + 1.0; 
    392398        zp2 = zz + 2.0; 
    … …  
    397403        zp7 = zz + 7.0; 
    398404        // 
    399         aa = (zz+1)/x/ravg; 
    400  
     405         
     406        //small QR limit - use Guinier approx 
     407        zp8 = zz+8.0; 
     408        if(x*ravg < 0.1) { 
     409                i_zero = scale*delrho*delrho*1.e8*4.*pi/3.*pow(ravg,3); 
     410                i_zero *= zp6*zp5*zp4/zp1/zp1/zp1;              //6th moment / 3rd moment 
     411                Rg2 = 3.*zp8*zp7/5./(zp1*zp1)*ravg*ravg; 
     412                pq = i_zero*exp(-x*x*Rg2/3.); 
     413                //pq += bkg;            //unlike the Igor code, the backgorund is added in the wrapper (above) 
     414                return(pq); 
     415        } 
     416        // 
     417 
     418        aa = (zz+1.0)/x/ravg; 
     419         
    401420        at1 = atan(1.0/aa); 
    402421        at2 = atan(2.0/aa); 
    … …  
    404423        //  calculations are performed to avoid  large # errors 
    405424        // - trick is to propogate the a^(z+7) term through the g1 
    406         // 
     425        //  
    407426        t1 = zp7*log10(aa) - zp1/2.0*log10(aa*aa+4.0); 
    408427        t2 = zp7*log10(aa) - zp3/2.0*log10(aa*aa+4.0); 
    … …  
    416435        v3 = -2.0*zp1*rt3*sin(zp2*at2); 
    417436        g1 = (v1+v2+v3); 
    418  
     437         
    419438        pq = log10(g1) - 6.0*log10(zp1) + 6.0*log10(ravg); 
    420         pq = pow(10,pq)*8*pi*pi*delrho*delrho; 
    421  
     439        pq = pow(10,pq)*8.0*pi*pi*delrho*delrho; 
     440         
    422441        // 
    423         // beta factor is not used here, but could be for the 
     442        // beta factor is not used here, but could be for the  
    424443        // decoupling approximation 
    425         // 
     444        //  
    426445        //      g11 = g1 
    427446        //      gd = -zp7*log(aa) 
    428447        //      g1 = log(g11) + gd 
    429         // 
     448        //                        
    430449        //      t1 = zp1*at1 
    431450        //      t2 = zp2*at1 
    432451        //      g2 = sin( t1 ) - zp1/sqrt(aa*aa+1)*cos( t2 ) 
    433452        //      g22 = g2*g2 
    434         //      beta = zp1*log(aa) - zp1*log(aa*aa+1) - g1 + log(g22) 
     453        //      beta = zp1*log(aa) - zp1*log(aa*aa+1) - g1 + log(g22)  
    435454        //      beta = 2*alog(beta) 
    436  
     455         
    437456        //re-normalize by the average volume 
    438457        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*ravg*ravg*ravg; 
    … …  
    440459        //scale, convert to cm^-1 
    441460        pq *= scale * 1.0e8; 
    442  
     461     
    443462    return(pq); 
    444463} 
    … …  
    451470        double pi,x; 
    452471        double scale,rad,pd,cont,bkg;           //my local names 
    453         double inten,h1,qw,qr,width,sig,averad3; 
    454  
     472        double inten,h1,qw,qr,width,sig,averad3,Rg2,slds,sld; 
     473         
    455474        pi = 4.0*atan(1.0); 
    456475        x= q; 
    457  
     476         
    458477        scale = dp[0]; 
    459478        rad = dp[1];            // radius (A) 
    460479        pd = dp[2];             //polydispersity of rectangular distribution 
    461         cont = dp[3];           // contrast (A^-2) 
    462         bkg = dp[4]; 
    463  
     480        slds = dp[3]; 
     481        sld = dp[4]; 
     482        cont = slds - sld;              // contrast (A^-2) 
     483        bkg = dp[5]; 
     484         
    464485        // as usual, poly = sig/ravg 
    465486        // for the rectangular distribution, sig = width/sqrt(3) 
    466487        // width is the HALF- WIDTH of the rectangular distrubution 
    467  
     488         
    468489        sig = pd*rad; 
    469490        width = sqrt(3.0)*sig; 
    470  
     491         
    471492        //x is the q-value 
    472493        qw = x*width; 
    473494        qr = x*rad; 
     495         
     496        // as for the numerical inatabilities at low QR, the function is calculating the sines and cosines 
     497        // just fine - the problem seems to be that the  
     498        // leading terms nearly cancel with the last term (the -6*qr... term), to within machine 
     499        // precision - the difference is on the order of 10^-20 
     500        // so just use the limiting Guiner value 
     501        if(qr<0.1) { 
     502                h1 = scale*cont*cont*1.e8*4.*pi/3.0*pow(rad,3); 
     503                h1 *=   (1. + 15.*pow(pd,2) + 27.*pow(pd,4) +27./7.*pow(pd,6) );                                //6th moment 
     504                h1 /= (1.+3.*pd*pd);    //3rd moment 
     505                Rg2 = 3.0/5.0*rad*rad*( 1.+28.*pow(pd,2)+126.*pow(pd,4)+108.*pow(pd,6)+27.*pow(pd,8) ); 
     506                Rg2 /= (1.+15.*pow(pd,2)+27.*pow(pd,4)+27./7.*pow(pd,6)); 
     507                h1 *= exp(-1./3.*Rg2*x*x); 
     508                h1 += bkg; 
     509                return(h1); 
     510        } 
     511         
     512        // normal calculation 
    474513        h1 = -0.5*qw + qr*qr*qw + (qw*qw*qw)/3.0; 
    475         h1 -= 5.0/2.0*cos(2*qr)*sin(qw)*cos(qw); 
    476         h1 += 0.5*qr*qr*cos(2*qr)*sin(2*qw); 
    477         h1 += 0.5*qw*qw*cos(2*qr)*sin(2*qw); 
    478         h1 += qw*qr*sin(2*qr)*cos(2*qw); 
     514        h1 -= 5.0/2.0*cos(2.0*qr)*sin(qw)*cos(qw); 
     515        h1 += 0.5*qr*qr*cos(2.0*qr)*sin(2.0*qw); 
     516        h1 += 0.5*qw*qw*cos(2.0*qr)*sin(2.0*qw); 
     517        h1 += qw*qr*sin(2.0*qr)*cos(2.0*qw); 
    479518        h1 += 3.0*qw*(cos(qr)*cos(qw))*(cos(qr)*cos(qw)); 
    480519        h1+= 3.0*qw*(sin(qr)*sin(qw))*(sin(qr)*sin(qw)); 
    481520        h1 -= 6.0*qr*cos(qr)*sin(qr)*cos(qw)*sin(qw); 
    482  
     521         
    483522        // calculate P(q) = <f^2> 
    484523        inten = 8.0*pi*pi*cont*cont/width/pow(x,7)*h1; 
    485  
     524         
    486525        // beta(q) would be calculated as 2/width/x/h1*h2*h2 
    487         // with 
     526        // with  
    488527        // h2 = 2*sin(x*rad)*sin(x*width)-x*rad*cos(x*rad)*sin(x*width)-x*width*sin(x*rad)*cos(x*width) 
    489  
     528         
    490529        // normalize to the average volume 
    491530        // <R^3> = ravg^3*(1+3*pd^2) 
    492531        // or... "zf"  = (1 + 3*p^2), which will be greater than one 
    493  
     532         
    494533        averad3 =  rad*rad*rad*(1.0+3.0*pd*pd); 
    495534        inten /= 4.0*pi/3.0*averad3; 
    … …  
    500539        // then add in the background 
    501540        inten += bkg; 
    502  
     541         
    503542        return(inten); 
    504543} 
    … …  
    514553        double va,vb,zi,yy,summ,inten; 
    515554        int nord=20,ii; 
    516  
     555         
    517556        pi = 4.0*atan(1.0); 
    518557        x= q; 
    519  
     558         
    520559        scale=dp[0]; 
    521560        rad=dp[1]; 
    … …  
    526565        delrho=rho-rhos; 
    527566        bkg=dp[5]; 
    528  
     567         
    529568        va = -4.0*sig + rad; 
    530         if (va<0) { 
    531                 va=0;           //to avoid numerical error when  va<0 (-ve q-value) 
     569        if (va<0.0) { 
     570                va=0.0;         //to avoid numerical error when  va<0 (-ve q-value) 
    532571        } 
    533572        vb = 4.0*sig +rad; 
    534  
     573         
    535574        summ = 0.0;             // initialize integral 
    536575        for(ii=0;ii<nord;ii+=1) { 
    … …  
    542581                yy *= yy; 
    543582                yy *= Gauss20Wt[ii] *  Gauss_distr(sig,rad,zi); 
    544  
     583                 
    545584                summ += yy;             //add to the running total of the quadrature 
    546585        } 
    547586        // calculate value of integral to return 
    548587        inten = (vb-va)/2.0*summ; 
    549  
     588         
    550589        //re-normalize by polydisperse sphere volume 
    551590        inten /= (4.0*pi/3.0*rad*rad*rad)*(1.0+3.0*pd*pd); 
    552  
     591         
    553592        inten *= 1.0e8; 
    554593        inten *= scale; 
    555594        inten += bkg; 
    556  
     595         
    557596    return(inten);      //scale, and add in the background 
    558597} 
    … …  
    567606        double va,vb,zi,yy,summ,inten; 
    568607        int nord=76,ii; 
    569  
     608         
    570609        pi = 4.0*atan(1.0); 
    571610        x= q; 
    572  
     611         
    573612        scale=dp[0]; 
    574613        rad=dp[1];              //rad is the median radius 
    … …  
    579618        delrho=rho-rhos; 
    580619        bkg=dp[5]; 
    581  
     620         
    582621        va = -3.5*sig + mu; 
    583622        va = exp(va); 
    584         if (va<0) { 
    585                 va=0;           //to avoid numerical error when  va<0 (-ve q-value) 
     623        if (va<0.0) { 
     624                va=0.0;         //to avoid numerical error when  va<0 (-ve q-value) 
    586625        } 
    587626        vb = 3.5*sig*(1.0+sig) +mu; 
    588627        vb = exp(vb); 
    589  
     628         
    590629        summ = 0.0;             // initialize integral 
    591630        for(ii=0;ii<nord;ii+=1) { 
    … …  
    597636                yy *= yy; 
    598637                yy *= Gauss76Wt[ii] *  LogNormal_distr(sig,mu,zi); 
    599  
     638                 
    600639                summ += yy;             //add to the running total of the quadrature 
    601640        } 
    602641        // calculate value of integral to return 
    603642        inten = (vb-va)/2.0*summ; 
    604  
     643         
    605644        //re-normalize by polydisperse sphere volume 
    606645        r3 = exp(3.0*mu + 9.0/2.0*sig*sig);             // <R^3> directly 
    607646        inten /= (4.0*pi/3.0*r3);               //polydisperse volume 
    608  
     647         
    609648        inten *= 1.0e8; 
    610649        inten *= scale; 
    611650        inten += bkg; 
    612  
     651         
    613652        return(inten); 
    614653} 
    … …  
    616655static double 
    617656LogNormal_distr(double sig, double mu, double pt) 
    618 { 
     657{        
    619658        double retval,pi; 
    620  
    621         pi = 4.0*atan(1.0); 
    622         retval = (1/ (sig*pt*sqrt(2.0*pi)) )*exp( -0.5*(log(pt) - mu)*(log(pt) - mu)/sig/sig ); 
     659         
     660        pi = 4.0*atan(1.0); 
     661        retval = (1.0/ (sig*pt*sqrt(2.0*pi)) )*exp( -0.5*(log(pt) - mu)*(log(pt) - mu)/sig/sig ); 
    623662        return(retval); 
    624663} 
    … …  
    626665static double 
    627666Gauss_distr(double sig, double avg, double pt) 
    628 { 
     667{        
    629668        double retval,Pi; 
    630  
     669         
    631670        Pi = 4.0*atan(1.0); 
    632671        retval = (1.0/ (sig*sqrt(2.0*Pi)) )*exp(-(avg-pt)*(avg-pt)/sig/sig/2.0); 
    … …  
    644683        double sld1,sld2,sld3,zp1,zp2,zp3,vpoly; 
    645684        double pi43,c1,c2,form,volume,arg1,arg2; 
    646  
     685         
    647686        pi = 4.0*atan(1.0); 
    648687        x= q; 
    649  
     688         
    650689        scale = dp[0]; 
    651690        corrad = dp[1]; 
    … …  
    656695        sld3 = dp[6]; 
    657696        bkg = dp[7]; 
    658  
    659         zz = (1.0/sig)*(1.0/sig) - 1.0; 
     697         
     698        zz = (1.0/sig)*(1.0/sig) - 1.0;    
    660699        shlrad = corrad + thick; 
    661700        drho1 = sld1-sld2;              //core-shell 
    … …  
    665704        zp3 = zz + 3.; 
    666705        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((corrad+thick),3); 
    667  
     706         
    668707        // the beta factor is not calculated 
    669708        // the calculated form factor <f^2> has units [length^2] 
    670709        // and must be multiplied by number density [l^-3] and the correct unit 
    671710        // conversion to get to absolute scale 
    672  
     711         
    673712        pi43=4.0/3.0*pi; 
    674713        pp=corrad/shlrad; 
    … …  
    676715        c1=drho1*volume; 
    677716        c2=drho2*volume; 
    678  
     717         
    679718        arg1 = x*shlrad*pp; 
    680719        arg2 = x*shlrad; 
    681  
     720         
    682721        form=pow(pp,6)*c1*c1*fnt2(arg1,zz); 
    683722        form += c2*c2*fnt2(arg2,zz); 
    684723        form += 2.0*c1*c2*fnt3(arg2,pp,zz); 
    685  
     724         
    686725        //convert the result to [cm^-1] 
    687  
     726         
    688727        //scale the result 
    689728        // - divide by the polydisperse volume, mult by 10^8 
    … …  
    691730        form *= 1.0e8; 
    692731        form *= scale; 
    693  
     732         
    694733        //add in the background 
    695734        form += bkg; 
    696  
     735         
    697736        return(form); 
    698737} 
    … …  
    706745{ 
    707746        double z1,z2,z3,u,ww,term1,term2,term3,ans; 
    708  
     747         
    709748        z1=zz+1.0; 
    710749        z2=zz+2.0; 
    … …  
    716755        term3=1.0+cos(z3*ww)/pow((1.0+4.0*u*u),(z3/2.0)); 
    717756        ans=(4.50/z1/pow(yy,6))*(z1*(1.0-term1-term2)+yy*yy*z2*term3); 
    718  
     757         
    719758        return(ans); 
    720759} 
    … …  
    726765double 
    727766fnt3(double yy, double pp, double zz) 
    728 { 
     767{        
    729768        double z1,z2,z3,yp,yn,up,un,vp,vn,term1,term2,term3,term4,term5,term6,ans; 
    730  
     769         
    731770        z1=zz+1.0; 
    732771        z2=zz+2.0; 
    … …  
    746785        ans=4.5/z1/pow(yy,6); 
    747786        ans *=(z1*(term1-term2)+yy*yy*pp*z2*(term3+term4)+z1*(term5-term6)); 
    748  
     787         
    749788        return(ans); 
    750789} 
    … …  
    771810        double v1,v2,n1,n2,qr1,qr2,b1,b2,sc1,sc2; 
    772811        int err; 
    773  
     812         
    774813        pi = 4.0*atan(1.0); 
    775814        x= q; 
    … …  
    782821        rhos = dp[6]; 
    783822        bgd = dp[7]; 
    784  
    785  
     823         
     824         
    786825        phi = phi1 + phi2; 
    787826        aa = r1/r2; 
    … …  
    792831        // calculate the PSF's here 
    793832        err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12); 
    794  
     833         
    795834        // /* do form factor calculations  */ 
    796  
     835         
    797836        v1 = 4.0*pi/3.0*r1*r1*r1; 
    798837        v2 = 4.0*pi/3.0*r2*r2*r2; 
    799  
     838         
    800839        n1 = phi1/v1; 
    801840        n2 = phi2/v2; 
    802  
     841         
    803842        qr1 = r1*x; 
    804843        qr2 = r2*x; 
    … …  
    821860        ///* convert I(1/A) to (1/cm)  */ 
    822861        inten *= 1.0e8; 
    823  
     862         
    824863        inten += bgd; 
    825  
     864         
    826865        return(inten); 
    827866} 
    … …  
    835874        double phi1,phi2,phr,a3; 
    836875        int err; 
    837  
     876         
    838877        pi = 4.0*atan(1.0); 
    839878        x= q; 
    … …  
    854893        // calculate the PSF's here 
    855894        err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12); 
    856  
     895         
    857896    return(psf11);      //scale, and add in the background 
    858897} 
    … …  
    866905        double phi1,phi2,phr,a3; 
    867906        int err; 
    868  
     907         
    869908        pi = 4.0*atan(1.0); 
    870909        x= q; 
    … …  
    885924        // calculate the PSF's here 
    886925        err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12); 
    887  
     926         
    888927    return(psf12);      //scale, and add in the background 
    889928} 
    … …  
    897936        double phi1,phi2,phr,a3; 
    898937        int err; 
    899  
     938         
    900939        pi = 4.0*atan(1.0); 
    901940        x= q; 
    902  
     941         
    903942        r2 = dp[0]; 
    904943        r1 = dp[1]; 
    … …  
    917956        // calculate the PSF's here 
    918957        err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12); 
    919  
     958         
    920959    return(psf22);      //scale, and add in the background 
    921960} 
    … …  
    925964{ 
    926965        //      variable qval,r2,nf2,aa,phi,&s11,&s22,&s12 
    927  
     966         
    928967        //   calculate constant terms 
    929968        double s1,s2,v,a3,v1,v2,g11,g12,g22,wmv,wmv3,wmv4; 
    930969        double a1,a2i,a2,b1,b2,b12,gm1,gm12; 
    931         double err=0,yy,ay,ay2,ay3,t1,t2,t3,f11,y2,y3,tt1,tt2,tt3; 
     970        double err=0.0,yy,ay,ay2,ay3,t1,t2,t3,f11,y2,y3,tt1,tt2,tt3; 
    932971        double c11,c22,c12,f12,f22,ttt1,ttt2,ttt3,ttt4,yl,y13; 
    933972        double t21,t22,t23,t31,t32,t33,t41,t42,yl3,wma3,y1; 
    934  
     973         
    935974        s2 = 2.0*r2; 
    936975        s1 = aa*s2; 
    … …  
    953992        gm1=(v1*a1+a3*v2*a2)*.5; 
    954993        gm12=2.*gm1*(1.-aa)/aa; 
    955         //c 
     994        //c   
    956995        //c   calculate the direct correlation functions and print results 
    957996        //c 
    958997        //      do 20 j=1,npts 
    959  
     998         
    960999        yy=qval*s2; 
    9611000        //c   calculate direct correlation functions 
    … …  
    9681007        t3=gm1*((4.*ay*ay2-24.*ay)*sin(ay)-(ay2*ay2-12.*ay2+24.)*cos(ay)+24.)/ay3; 
    9691008        f11=24.*v1*(t1+t2+t3)/ay3; 
    970  
     1009         
    9711010        //c ------c22 
    9721011        y2=yy*yy; 
    … …  
    9761015        tt3=gm1*((4.*y3-24.*yy)*sin(yy)-(y2*y2-12.*y2+24.)*cos(yy)+24.)/ay3; 
    9771016        f22=24.*v2*(tt1+tt2+tt3)/y3; 
    978  
     1017         
    9791018        //c   -----c12 
    9801019        yl=.5*yy*(1.-aa); 
    … …  
    9961035        ttt4=a1*(t41+t42)/y1; 
    9971036        f12=ttt1+24.*v*sqrt(nf2)*sqrt(1.-nf2)*a3*(ttt2+ttt3+ttt4)/(nf2+(1.-nf2)*a3); 
    998  
     1037         
    9991038        c11=f11; 
    10001039        c22=f22; 
    10011040        c12=f12; 
    10021041        *s11=1./(1.+c11-(c12)*c12/(1.+c22)); 
    1003         *s22=1./(1.+c22-(c12)*c12/(1.+c11)); 
    1004         *s12=-c12/((1.+c11)*(1.+c22)-(c12)*(c12)); 
    1005  
     1042        *s22=1./(1.+c22-(c12)*c12/(1.+c11));  
     1043        *s12=-c12/((1.+c11)*(1.+c22)-(c12)*(c12));    
     1044         
    10061045        return(err); 
    10071046} 
    … …  
    10161055 // bragg peaks arise naturally from the periodicity of the sample 
    10171056 // resolution smeared version gives he most appropriate view of the model 
    1018  
     1057  
    10191058        Warning: 
    10201059 The call to WaveData() below returns a pointer to the middle 
    … …  
    10311070        double fval,voli,ri,sldi; 
    10321071        double pi; 
    1033  
    1034         pi = 4.0*atan(1.0); 
    1035  
     1072         
     1073        pi = 4.0*atan(1.0); 
     1074         
    10361075        x= q; 
    10371076        scale = dp[0]; 
    … …  
    10431082        num = dp[6]; 
    10441083        bkg = dp[7]; 
    1045  
     1084         
    10461085        //calculate with a loop, two shells at a time 
    1047  
     1086         
    10481087        ii=0; 
    1049         fval=0; 
    1050  
     1088        fval=0.0; 
     1089         
    10511090        do { 
    10521091                ri = rcore + (double)ii*ts + (double)ii*tw; 
    1053                 voli = 4*pi/3*ri*ri*ri; 
     1092                voli = 4.0*pi/3.0*ri*ri*ri; 
    10541093                sldi = rhocore-rhoshel; 
    10551094                fval += voli*sldi*F_func(ri*x); 
    10561095                ri += ts; 
    1057                 voli = 4*pi/3*ri*ri*ri; 
     1096                voli = 4.0*pi/3.0*ri*ri*ri; 
    10581097                sldi = rhoshel-rhocore; 
    10591098                fval += voli*sldi*F_func(ri*x); 
    10601099                ii+=1;          //do 2 layers at a time 
    10611100        } while(ii<=num-1);  //change to make 0 < num < 2 correspond to unilamellar vesicles (C. Glinka, 11/24/03) 
    1062  
     1101         
    10631102        fval *= fval;           //square it 
    10641103        fval /= voli;           //normalize by the overall volume 
    1065         fval *= scale*1e8; 
     1104        fval *= scale*1.0e8; 
    10661105        fval += bkg; 
    1067  
     1106         
    10681107        return(fval); 
    10691108} 
    … …  
    10951134        double fval,ri,pi; 
    10961135        double avg,pd,zz,lo,hi,r1,r2,d1,d2,distr; 
    1097  
    1098         pi = 4.0*atan(1.0); 
     1136         
     1137        pi = 4.0*atan(1.0);      
    10991138        x= q; 
    1100  
     1139         
    11011140        scale = dp[0]; 
    11021141        avg = dp[1];            // average (total) outer radius 
    … …  
    11081147        rhoshel = dp[7]; 
    11091148        bkg = dp[8]; 
    1110  
     1149         
    11111150        zz = (1.0/pd)*(1.0/pd)-1.0; 
    1112  
     1151         
    11131152        //max radius set to be 5 std deviations past mean 
    11141153        hi = avg + pd*avg*5.0; 
    11151154        lo = avg - pd*avg*5.0; 
    1116  
     1155         
    11171156        maxPairs = trunc( (hi-rcore+tw)/(ts+tw) ); 
    11181157        minPairs = trunc( (lo-rcore+tw)/(ts+tw) ); 
    11191158        minPairs = (minPairs < 1) ? 1 : minPairs;       // need a minimum of one 
    1120  
     1159         
    11211160        ii=minPairs; 
    1122         fval=0; 
    1123         d1 = 0; 
    1124         d2 = 0; 
    1125         r1 = 0; 
    1126         r2 = 0; 
    1127         distr = 0; 
    1128         first = 1; 
     1161        fval=0.0; 
     1162        d1 = 0.0; 
     1163        d2 = 0.0; 
     1164        r1 = 0.0; 
     1165        r2 = 0.0; 
     1166        distr = 0.0; 
     1167        first = 1.0; 
    11291168        do { 
    11301169                //make the current values old 
    11311170                r1 = r2; 
    11321171                d1 = d2; 
    1133  
     1172                 
    11341173                ri = (double)ii*(ts+tw) - tw + rcore; 
    11351174                fval += SchulzPoint(ri,avg,zz) * MultiShellGuts(x, rcore, ts, tw, rhocore, rhoshel, ii) * (4*pi/3*ri*ri*ri); 
    … …  
    11431182                first = 0; 
    11441183        } while(ii<=maxPairs); 
    1145  
    1146         fval /= 4*pi/3*avg*avg*avg;             //normalize by the overall volume 
     1184         
     1185        fval /= 4.0*pi/3.0*avg*avg*avg;         //normalize by the overall volume 
    11471186        fval /= distr; 
    11481187        fval *= scale; 
    11491188        fval += bkg; 
    1150  
     1189         
    11511190        return(fval); 
    11521191} 
    … …  
    11541193double 
    11551194MultiShellGuts(double x,double rcore,double ts,double tw,double rhocore,double rhoshel,int num) { 
    1156  
     1195         
    11571196    double ri,sldi,fval,voli,pi; 
    11581197    int ii; 
    1159  
     1198     
    11601199        pi = 4.0*atan(1.0); 
    11611200    ii=0; 
    1162     fval=0; 
    1163  
     1201    fval=0.0; 
     1202     
    11641203    do { 
    11651204        ri = rcore + (double)ii*ts + (double)ii*tw; 
    1166         voli = 4*pi/3*ri*ri*ri; 
     1205        voli = 4.0*pi/3.0*ri*ri*ri; 
    11671206        sldi = rhocore-rhoshel; 
    11681207        fval += voli*sldi*F_func(ri*x); 
    11691208        ri += ts; 
    1170         voli = 4*pi/3*ri*ri*ri; 
     1209        voli = 4.0*pi/3.0*ri*ri*ri; 
    11711210        sldi = rhoshel-rhocore; 
    11721211        fval += voli*sldi*F_func(ri*x); 
    11731212        ii+=1;          //do 2 layers at a time 
    11741213    } while(ii<=num-1);  //change to make 0 < num < 2 correspond to unilamellar vesicles (C. Glinka, 11/24/03) 
    1175  
     1214     
    11761215    fval *= fval; 
    11771216    fval /= voli; 
    1178     fval *= 1e8; 
    1179  
     1217    fval *= 1.0e8; 
     1218     
    11801219    return(fval);       // this result still needs to be multiplied by scale and have background added 
    11811220} 
    … …  
    11831222static double 
    11841223SchulzPoint(double x, double avg, double zz) { 
    1185  
     1224         
    11861225    double dr; 
    1187  
    1188     dr = zz*log(x) - gammln(zz+1)+(zz+1)*log((zz+1)/avg)-(x/avg*(zz+1)); 
     1226     
     1227    dr = zz*log(x) - gammln(zz+1.0)+(zz+1.0)*log((zz+1.0)/avg)-(x/avg*(zz+1.0)); 
    11891228    return (exp(dr)); 
    11901229} 
    … …  
    11921231static double 
    11931232gammln(double xx) { 
    1194  
     1233         
    11951234    double x,y,tmp,ser; 
    11961235    static double cof[6]={76.18009172947146,-86.50532032941677, 
    … …  
    11981237                0.1208650973866179e-2,-0.5395239384953e-5}; 
    11991238    int j; 
    1200  
     1239         
    12011240    y=x=xx; 
    12021241    tmp=x+5.5; 
    … …  
    12101249F_func(double qr) { 
    12111250        double sc; 
    1212         if (qr == 0){ 
     1251        if (qr == 0.0){ 
    12131252                sc = 1.0; 
    12141253        }else{ 
    1215                 sc=(3*(sin(qr) - qr*cos(qr))/(qr*qr*qr)); 
     1254                sc=(3.0*(sin(qr) - qr*cos(qr))/(qr*qr*qr)); 
    12161255        } 
    12171256        return sc; 
    12181257} 
    12191258 
     1259double 
     1260OneShell(double dp[], double q) 
     1261{ 
     1262        // variables are: 
     1263        //[0] scale factor 
     1264        //[1] radius of core [] 
     1265        //[2] SLD of the core   [-2] 
     1266        //[3] thickness of the shell    [] 
     1267        //[4] SLD of the shell 
     1268        //[5] SLD of the solvent 
     1269        //[6] background        [cm-1] 
     1270         
     1271        double x,pi; 
     1272        double scale,rcore,thick,rhocore,rhoshel,rhosolv,bkg;           //my local names 
     1273        double bes,f,vol,qr,contr,f2; 
     1274         
     1275        pi = 4.0*atan(1.0); 
     1276        x=q; 
     1277         
     1278        scale = dp[0]; 
     1279        rcore = dp[1]; 
     1280        rhocore = dp[2]; 
     1281        thick = dp[3]; 
     1282        rhoshel = dp[4]; 
     1283        rhosolv = dp[5]; 
     1284        bkg = dp[6]; 
     1285         
     1286        // core first, then add in shell 
     1287        qr=x*rcore; 
     1288        contr = rhocore-rhoshel; 
     1289        if(qr == 0){ 
     1290                bes = 1.0; 
     1291        }else{ 
     1292                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1293        } 
     1294        vol = 4.0*pi/3.0*rcore*rcore*rcore; 
     1295        f = vol*bes*contr; 
     1296        //now the shell 
     1297        qr=x*(rcore+thick); 
     1298        contr = rhoshel-rhosolv; 
     1299        if(qr == 0){ 
     1300                bes = 1.0; 
     1301        }else{ 
     1302                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1303        } 
     1304        vol = 4.0*pi/3.0*pow((rcore+thick),3); 
     1305        f += vol*bes*contr; 
     1306         
     1307        // normalize to particle volume and rescale from [-1] to [cm-1] 
     1308        f2 = f*f/vol*1.0e8; 
     1309         
     1310        //scale if desired 
     1311        f2 *= scale; 
     1312        // then add in the background 
     1313        f2 += bkg; 
     1314         
     1315        return(f2); 
     1316} 
     1317 
     1318double 
     1319TwoShell(double dp[], double q) 
     1320{ 
     1321        // variables are: 
     1322        //[0] scale factor 
     1323        //[1] radius of core [] 
     1324        //[2] SLD of the core   [-2] 
     1325        //[3] thickness of shell 1 [] 
     1326        //[4] SLD of shell 1 
     1327        //[5] thickness of shell 2 [] 
     1328        //[6] SLD of shell 2 
     1329        //[7] SLD of the solvent 
     1330        //[8] background        [cm-1] 
     1331         
     1332        double x,pi; 
     1333        double scale,rcore,thick1,rhocore,rhoshel1,rhosolv,bkg;         //my local names 
     1334        double bes,f,vol,qr,contr,f2; 
     1335        double rhoshel2,thick2; 
     1336         
     1337        pi = 4.0*atan(1.0); 
     1338        x=q; 
     1339         
     1340        scale = dp[0]; 
     1341        rcore = dp[1]; 
     1342        rhocore = dp[2]; 
     1343        thick1 = dp[3]; 
     1344        rhoshel1 = dp[4]; 
     1345        thick2 = dp[5]; 
     1346        rhoshel2 = dp[6];        
     1347        rhosolv = dp[7]; 
     1348        bkg = dp[8]; 
     1349                // core first, then add in shells 
     1350        qr=x*rcore; 
     1351        contr = rhocore-rhoshel1; 
     1352        if(qr == 0){ 
     1353                bes = 1.0; 
     1354        }else{ 
     1355                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1356        } 
     1357        vol = 4.0*pi/3.0*rcore*rcore*rcore; 
     1358        f = vol*bes*contr; 
     1359        //now the shell (1) 
     1360        qr=x*(rcore+thick1); 
     1361        contr = rhoshel1-rhoshel2; 
     1362        if(qr == 0){ 
     1363                bes = 1.0; 
     1364        }else{ 
     1365                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1366        } 
     1367        vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1); 
     1368        f += vol*bes*contr; 
     1369        //now the shell (2) 
     1370        qr=x*(rcore+thick1+thick2); 
     1371        contr = rhoshel2-rhosolv; 
     1372        if(qr == 0){ 
     1373                bes = 1.0; 
     1374        }else{ 
     1375                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1376        } 
     1377        vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2); 
     1378        f += vol*bes*contr; 
     1379         
     1380                 
     1381        // normalize to particle volume and rescale from [-1] to [cm-1] 
     1382        f2 = f*f/vol*1.0e8; 
     1383         
     1384        //scale if desired 
     1385        f2 *= scale; 
     1386        // then add in the background 
     1387        f2 += bkg; 
     1388         
     1389        return(f2); 
     1390} 
     1391 
     1392double 
     1393ThreeShell(double dp[], double q) 
     1394{ 
     1395        // variables are: 
     1396        //[0] scale factor 
     1397        //[1] radius of core [] 
     1398        //[2] SLD of the core   [-2] 
     1399        //[3] thickness of shell 1 [] 
     1400        //[4] SLD of shell 1 
     1401        //[5] thickness of shell 2 [] 
     1402        //[6] SLD of shell 2 
     1403        //[7] thickness of shell 3 
     1404        //[8] SLD of shell 3 
     1405        //[9] SLD of solvent 
     1406        //[10] background       [cm-1] 
     1407         
     1408        double x,pi; 
     1409        double scale,rcore,thick1,rhocore,rhoshel1,rhosolv,bkg;         //my local names 
     1410        double bes,f,vol,qr,contr,f2; 
     1411        double rhoshel2,thick2,rhoshel3,thick3; 
     1412         
     1413        pi = 4.0*atan(1.0); 
     1414        x=q; 
     1415         
     1416        scale = dp[0]; 
     1417        rcore = dp[1]; 
     1418        rhocore = dp[2]; 
     1419        thick1 = dp[3]; 
     1420        rhoshel1 = dp[4]; 
     1421        thick2 = dp[5]; 
     1422        rhoshel2 = dp[6];        
     1423        thick3 = dp[7]; 
     1424        rhoshel3 = dp[8];        
     1425        rhosolv = dp[9]; 
     1426        bkg = dp[10]; 
     1427         
     1428                // core first, then add in shells 
     1429        qr=x*rcore; 
     1430        contr = rhocore-rhoshel1; 
     1431        if(qr == 0){ 
     1432                bes = 1.0; 
     1433        }else{ 
     1434                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1435        } 
     1436        vol = 4.0*pi/3.0*rcore*rcore*rcore; 
     1437        f = vol*bes*contr; 
     1438        //now the shell (1) 
     1439        qr=x*(rcore+thick1); 
     1440        contr = rhoshel1-rhoshel2; 
     1441        if(qr == 0){ 
     1442                bes = 1.0; 
     1443        }else{ 
     1444                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1445        } 
     1446        vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1); 
     1447        f += vol*bes*contr; 
     1448        //now the shell (2) 
     1449        qr=x*(rcore+thick1+thick2); 
     1450        contr = rhoshel2-rhoshel3; 
     1451        if(qr == 0){ 
     1452                bes = 1.0; 
     1453        }else{ 
     1454                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1455        } 
     1456        vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2); 
     1457        f += vol*bes*contr; 
     1458        //now the shell (3) 
     1459        qr=x*(rcore+thick1+thick2+thick3); 
     1460        contr = rhoshel3-rhosolv; 
     1461        if(qr == 0){ 
     1462                bes = 1.0; 
     1463        }else{ 
     1464                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1465        } 
     1466        vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3); 
     1467        f += vol*bes*contr; 
     1468                 
     1469        // normalize to particle volume and rescale from [-1] to [cm-1] 
     1470        f2 = f*f/vol*1.0e8; 
     1471         
     1472        //scale if desired 
     1473        f2 *= scale; 
     1474        // then add in the background 
     1475        f2 += bkg; 
     1476         
     1477        return(f2); 
     1478} 
     1479 
     1480double 
     1481FourShell(double dp[], double q) 
     1482{ 
     1483        // variables are: 
     1484        //[0] scale factor 
     1485        //[1] radius of core [] 
     1486        //[2] SLD of the core   [-2] 
     1487        //[3] thickness of shell 1 [] 
     1488        //[4] SLD of shell 1 
     1489        //[5] thickness of shell 2 [] 
     1490        //[6] SLD of shell 2 
     1491        //[7] thickness of shell 3 
     1492        //[8] SLD of shell 3 
     1493        //[9] thickness of shell 3 
     1494        //[10] SLD of shell 3 
     1495        //[11] SLD of solvent 
     1496        //[12] background       [cm-1] 
     1497         
     1498        double x,pi; 
     1499        double scale,rcore,thick1,rhocore,rhoshel1,rhosolv,bkg;         //my local names 
     1500        double bes,f,vol,qr,contr,f2; 
     1501        double rhoshel2,thick2,rhoshel3,thick3,rhoshel4,thick4; 
     1502         
     1503        pi = 4.0*atan(1.0); 
     1504        x=q; 
     1505         
     1506        scale = dp[0]; 
     1507        rcore = dp[1]; 
     1508        rhocore = dp[2]; 
     1509        thick1 = dp[3]; 
     1510        rhoshel1 = dp[4]; 
     1511        thick2 = dp[5]; 
     1512        rhoshel2 = dp[6];        
     1513        thick3 = dp[7]; 
     1514        rhoshel3 = dp[8]; 
     1515        thick4 = dp[9]; 
     1516        rhoshel4 = dp[10];       
     1517        rhosolv = dp[11]; 
     1518        bkg = dp[12]; 
     1519         
     1520                // core first, then add in shells 
     1521        qr=x*rcore; 
     1522        contr = rhocore-rhoshel1; 
     1523        if(qr == 0){ 
     1524                bes = 1.0; 
     1525        }else{ 
     1526                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1527        } 
     1528        vol = 4.0*pi/3.0*rcore*rcore*rcore; 
     1529        f = vol*bes*contr; 
     1530        //now the shell (1) 
     1531        qr=x*(rcore+thick1); 
     1532        contr = rhoshel1-rhoshel2; 
     1533        if(qr == 0){ 
     1534                bes = 1.0; 
     1535        }else{ 
     1536                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1537        } 
     1538        vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1); 
     1539        f += vol*bes*contr; 
     1540        //now the shell (2) 
     1541        qr=x*(rcore+thick1+thick2); 
     1542        contr = rhoshel2-rhoshel3; 
     1543        if(qr == 0){ 
     1544                bes = 1.0; 
     1545        }else{ 
     1546                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1547        } 
     1548        vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2); 
     1549        f += vol*bes*contr; 
     1550        //now the shell (3) 
     1551        qr=x*(rcore+thick1+thick2+thick3); 
     1552        contr = rhoshel3-rhoshel4; 
     1553        if(qr == 0){ 
     1554                bes = 1.0; 
     1555        }else{ 
     1556                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1557        } 
     1558        vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3); 
     1559        f += vol*bes*contr; 
     1560        //now the shell (4) 
     1561        qr=x*(rcore+thick1+thick2+thick3+thick4); 
     1562        contr = rhoshel4-rhosolv; 
     1563        if(qr == 0){ 
     1564                bes = 1.0; 
     1565        }else{ 
     1566                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1567        } 
     1568        vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4); 
     1569        f += vol*bes*contr; 
     1570         
     1571                 
     1572        // normalize to particle volume and rescale from [-1] to [cm-1] 
     1573        f2 = f*f/vol*1.0e8; 
     1574         
     1575        //scale if desired 
     1576        f2 *= scale; 
     1577        // then add in the background 
     1578        f2 += bkg; 
     1579         
     1580        return(f2); 
     1581} 
     1582 
     1583double 
     1584PolyOneShell(double dp[], double x) 
     1585{ 
     1586        double scale,rcore,thick,rhocore,rhoshel,rhosolv,bkg,pd,zz;             //my local names 
     1587        double va,vb,summ,yyy,zi; 
     1588        double answer,zp1,zp2,zp3,vpoly,range,temp_1sf[7],pi; 
     1589        int nord=76,ii; 
     1590         
     1591        pi = 4.0*atan(1.0); 
     1592         
     1593        scale = dp[0]; 
     1594        rcore = dp[1]; 
     1595        pd = dp[2]; 
     1596        rhocore = dp[3]; 
     1597        thick = dp[4]; 
     1598        rhoshel = dp[5]; 
     1599        rhosolv = dp[6]; 
     1600        bkg = dp[7]; 
     1601                 
     1602        zz = (1.0/pd)*(1.0/pd)-1.0;             //polydispersity of the core only 
     1603         
     1604        range = 8.0;                    //std deviations for the integration 
     1605        va = rcore*(1.0-range*pd); 
     1606        if (va<0.0) { 
     1607                va=0.0;         //otherwise numerical error when pd >= 0.3, making a<0 
     1608        } 
     1609        if (pd>0.3) { 
     1610                range = range + (pd-0.3)*18.0;          //stretch upper range to account for skewed tail 
     1611        } 
     1612        vb = rcore*(1.0+range*pd);              // is this far enough past avg radius? 
     1613         
     1614        //temp set scale=1 and bkg=0 for quadrature calc 
     1615        temp_1sf[0] = 1.0; 
     1616        temp_1sf[1] = dp[1];    //the core radius will be changed in the loop 
     1617        temp_1sf[2] = dp[3]; 
     1618        temp_1sf[3] = dp[4]; 
     1619        temp_1sf[4] = dp[5]; 
     1620        temp_1sf[5] = dp[6]; 
     1621        temp_1sf[6] = 0.0; 
     1622         
     1623        summ = 0.0;             // initialize integral 
     1624        for(ii=0;ii<nord;ii+=1) { 
     1625                // calculate Gauss points on integration interval (r-value for evaluation) 
     1626                zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0; 
     1627                temp_1sf[1] = zi; 
     1628                yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * OneShell(temp_1sf,x); 
     1629                //un-normalize by volume 
     1630                yyy *= 4.0*pi/3.0*pow((zi+thick),3); 
     1631                summ += yyy;            //add to the running total of the quadrature 
     1632        } 
     1633        // calculate value of integral to return 
     1634        answer = (vb-va)/2.0*summ; 
     1635         
     1636        //re-normalize by the average volume 
     1637        zp1 = zz + 1.0; 
     1638        zp2 = zz + 2.0; 
     1639        zp3 = zz + 3.0; 
     1640        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick),3); 
     1641        answer /= vpoly; 
     1642//scale 
     1643        answer *= scale; 
     1644// add in the background 
     1645        answer += bkg; 
     1646                 
     1647    return(answer); 
     1648} 
     1649 
     1650double 
     1651PolyTwoShell(double dp[], double x) 
     1652{ 
     1653        double scale,rcore,rhocore,rhosolv,bkg,pd,zz;           //my local names 
     1654        double va,vb,summ,yyy,zi; 
     1655        double answer,zp1,zp2,zp3,vpoly,range,temp_2sf[9],pi; 
     1656        int nord=76,ii; 
     1657        double thick1,thick2; 
     1658        double rhoshel1,rhoshel2; 
     1659         
     1660        scale = dp[0]; 
     1661        rcore = dp[1]; 
     1662        pd = dp[2]; 
     1663        rhocore = dp[3]; 
     1664        thick1 = dp[4]; 
     1665        rhoshel1 = dp[5]; 
     1666        thick2 = dp[6]; 
     1667        rhoshel2 = dp[7]; 
     1668        rhosolv = dp[8]; 
     1669        bkg = dp[9]; 
     1670         
     1671        pi = 4.0*atan(1.0); 
     1672                 
     1673        zz = (1.0/pd)*(1.0/pd)-1.0;             //polydispersity of the core only 
     1674         
     1675        range = 8.0;                    //std deviations for the integration 
     1676        va = rcore*(1.0-range*pd); 
     1677        if (va<0.0) { 
     1678                va=0.0;         //otherwise numerical error when pd >= 0.3, making a<0 
     1679        } 
     1680        if (pd>0.3) { 
     1681                range = range + (pd-0.3)*18.0;          //stretch upper range to account for skewed tail 
     1682        } 
     1683        vb = rcore*(1.0+range*pd);              // is this far enough past avg radius? 
     1684         
     1685        //temp set scale=1 and bkg=0 for quadrature calc 
     1686        temp_2sf[0] = 1.0; 
     1687        temp_2sf[1] = dp[1];            //the core radius will be changed in the loop 
     1688        temp_2sf[2] = dp[3]; 
     1689        temp_2sf[3] = dp[4]; 
     1690        temp_2sf[4] = dp[5]; 
     1691        temp_2sf[5] = dp[6]; 
     1692        temp_2sf[6] = dp[7]; 
     1693        temp_2sf[7] = dp[8]; 
     1694        temp_2sf[8] = 0.0; 
     1695         
     1696        summ = 0.0;             // initialize integral 
     1697        for(ii=0;ii<nord;ii+=1) { 
     1698                // calculate Gauss points on integration interval (r-value for evaluation) 
     1699                zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0; 
     1700                temp_2sf[1] = zi; 
     1701                yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * TwoShell(temp_2sf,x); 
     1702                //un-normalize by volume 
     1703                yyy *= 4.0*pi/3.0*pow((zi+thick1+thick2),3); 
     1704                summ += yyy;            //add to the running total of the quadrature 
     1705        } 
     1706        // calculate value of integral to return 
     1707        answer = (vb-va)/2.0*summ; 
     1708         
     1709        //re-normalize by the average volume 
     1710        zp1 = zz + 1.0; 
     1711        zp2 = zz + 2.0; 
     1712        zp3 = zz + 3.0; 
     1713        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick1+thick2),3); 
     1714        answer /= vpoly; 
     1715//scale 
     1716        answer *= scale; 
     1717// add in the background 
     1718        answer += bkg; 
     1719                 
     1720    return(answer); 
     1721} 
     1722 
     1723double 
     1724PolyThreeShell(double dp[], double x) 
     1725{ 
     1726        double scale,rcore,rhocore,rhosolv,bkg,pd,zz;           //my local names 
     1727        double va,vb,summ,yyy,zi; 
     1728        double answer,zp1,zp2,zp3,vpoly,range,temp_3sf[11],pi; 
     1729        int nord=76,ii; 
     1730        double thick1,thick2,thick3; 
     1731        double rhoshel1,rhoshel2,rhoshel3; 
     1732         
     1733        scale = dp[0]; 
     1734        rcore = dp[1]; 
     1735        pd = dp[2]; 
     1736        rhocore = dp[3]; 
     1737        thick1 = dp[4]; 
     1738        rhoshel1 = dp[5]; 
     1739        thick2 = dp[6]; 
     1740        rhoshel2 = dp[7]; 
     1741        thick3 = dp[8]; 
     1742        rhoshel3 = dp[9]; 
     1743        rhosolv = dp[10]; 
     1744        bkg = dp[11]; 
     1745         
     1746        pi = 4.0*atan(1.0); 
     1747                 
     1748        zz = (1.0/pd)*(1.0/pd)-1.0;             //polydispersity of the core only 
     1749         
     1750        range = 8.0;                    //std deviations for the integration 
     1751        va = rcore*(1.0-range*pd); 
     1752        if (va<0) { 
     1753                va=0;           //otherwise numerical error when pd >= 0.3, making a<0 
     1754        } 
     1755        if (pd>0.3) { 
     1756                range = range + (pd-0.3)*18.0;          //stretch upper range to account for skewed tail 
     1757        } 
     1758        vb = rcore*(1.0+range*pd);              // is this far enough past avg radius? 
     1759         
     1760        //temp set scale=1 and bkg=0 for quadrature calc 
     1761        temp_3sf[0] = 1.0; 
     1762        temp_3sf[1] = dp[1];            //the core radius will be changed in the loop 
     1763        temp_3sf[2] = dp[3]; 
     1764        temp_3sf[3] = dp[4]; 
     1765        temp_3sf[4] = dp[5]; 
     1766        temp_3sf[5] = dp[6]; 
     1767        temp_3sf[6] = dp[7]; 
     1768        temp_3sf[7] = dp[8]; 
     1769        temp_3sf[8] = dp[9]; 
     1770        temp_3sf[9] = dp[10]; 
     1771        temp_3sf[10] = 0.0; 
     1772         
     1773        summ = 0.0;             // initialize integral 
     1774        for(ii=0;ii<nord;ii+=1) { 
     1775                // calculate Gauss points on integration interval (r-value for evaluation) 
     1776                zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0; 
     1777                temp_3sf[1] = zi; 
     1778                yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * ThreeShell(temp_3sf,x); 
     1779                //un-normalize by volume 
     1780                yyy *= 4.0*pi/3.0*pow((zi+thick1+thick2+thick3),3); 
     1781                summ += yyy;            //add to the running total of the quadrature 
     1782        } 
     1783        // calculate value of integral to return 
     1784        answer = (vb-va)/2.0*summ; 
     1785         
     1786        //re-normalize by the average volume 
     1787        zp1 = zz + 1.0; 
     1788        zp2 = zz + 2.0; 
     1789        zp3 = zz + 3.0; 
     1790        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick1+thick2+thick3),3); 
     1791        answer /= vpoly; 
     1792//scale 
     1793        answer *= scale; 
     1794// add in the background 
     1795        answer += bkg; 
     1796                 
     1797    return(answer); 
     1798} 
     1799 
     1800double 
     1801PolyFourShell(double dp[], double x) 
     1802{ 
     1803        double scale,rcore,rhocore,rhosolv,bkg,pd,zz;           //my local names 
     1804        double va,vb,summ,yyy,zi; 
     1805        double answer,zp1,zp2,zp3,vpoly,range,temp_4sf[13],pi; 
     1806        int nord=76,ii; 
     1807        double thick1,thick2,thick3,thick4; 
     1808        double rhoshel1,rhoshel2,rhoshel3,rhoshel4; 
     1809         
     1810        scale = dp[0]; 
     1811        rcore = dp[1]; 
     1812        pd = dp[2]; 
     1813        rhocore = dp[3]; 
     1814        thick1 = dp[4]; 
     1815        rhoshel1 = dp[5]; 
     1816        thick2 = dp[6]; 
     1817        rhoshel2 = dp[7]; 
     1818        thick3 = dp[8]; 
     1819        rhoshel3 = dp[9]; 
     1820        thick4 = dp[10]; 
     1821        rhoshel4 = dp[11]; 
     1822        rhosolv = dp[12]; 
     1823        bkg = dp[13]; 
     1824         
     1825        pi = 4.0*atan(1.0); 
     1826                 
     1827        zz = (1.0/pd)*(1.0/pd)-1.0;             //polydispersity of the core only 
     1828         
     1829        range = 8.0;                    //std deviations for the integration 
     1830        va = rcore*(1.0-range*pd); 
     1831        if (va<0) { 
     1832                va=0;           //otherwise numerical error when pd >= 0.3, making a<0 
     1833        } 
     1834        if (pd>0.3) { 
     1835                range = range + (pd-0.3)*18.0;          //stretch upper range to account for skewed tail 
     1836        } 
     1837        vb = rcore*(1.0+range*pd);              // is this far enough past avg radius? 
     1838         
     1839        //temp set scale=1 and bkg=0 for quadrature calc 
     1840        temp_4sf[0] = 1.0; 
     1841        temp_4sf[1] = dp[1];            //the core radius will be changed in the loop 
     1842        temp_4sf[2] = dp[3]; 
     1843        temp_4sf[3] = dp[4]; 
     1844        temp_4sf[4] = dp[5]; 
     1845        temp_4sf[5] = dp[6]; 
     1846        temp_4sf[6] = dp[7]; 
     1847        temp_4sf[7] = dp[8]; 
     1848        temp_4sf[8] = dp[9]; 
     1849        temp_4sf[9] = dp[10]; 
     1850        temp_4sf[10] = dp[11]; 
     1851        temp_4sf[11] = dp[12]; 
     1852        temp_4sf[12] = 0.0; 
     1853                 
     1854        summ = 0.0;             // initialize integral 
     1855        for(ii=0;ii<nord;ii+=1) { 
     1856                // calculate Gauss points on integration interval (r-value for evaluation) 
     1857                zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0; 
     1858                temp_4sf[1] = zi; 
     1859                yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * FourShell(temp_4sf,x); 
     1860                //un-normalize by volume 
     1861                yyy *= 4.0*pi/3.0*pow((zi+thick1+thick2+thick3+thick4),3); 
     1862                summ += yyy;            //add to the running total of the quadrature 
     1863        } 
     1864        // calculate value of integral to return 
     1865        answer = (vb-va)/2.0*summ; 
     1866         
     1867        //re-normalize by the average volume 
     1868        zp1 = zz + 1.0; 
     1869        zp2 = zz + 2.0; 
     1870        zp3 = zz + 3.0; 
     1871        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick1+thick2+thick3+thick4),3); 
     1872        answer /= vpoly; 
     1873//scale 
     1874        answer *= scale; 
     1875// add in the background 
     1876        answer += bkg; 
     1877                 
     1878    return(answer); 
     1879} 
     1880 
     1881 
     1882/*      BCC_ParaCrystal  :  calculates the form factor of a Triaxial Ellipsoid at the given x-value p->x 
     1883 
     1884Uses 150 pt Gaussian quadrature for both integrals 
     1885 
     1886*/ 
     1887double 
     1888BCC_ParaCrystal(double w[], double x) 
     1889{ 
     1890        int i,j; 
     1891        double Pi; 
     1892        double scale,Dnn,gg,Rad,contrast,background,latticeScale,sld,sldSolv;           //local variables of coefficient wave 
     1893        int nordi=150;                  //order of integration 
     1894        int nordj=150; 
     1895        double va,vb;           //upper and lower integration limits 
     1896        double summ,zi,yyy,answer;                      //running tally of integration 
     1897        double summj,vaj,vbj,zij;                       //for the inner integration 
     1898         
     1899        Pi = 4.0*atan(1.0); 
     1900        va = 0.0; 
     1901        vb = 2.0*Pi;            //orintational average, outer integral 
     1902        vaj = 0.0; 
     1903        vbj = Pi;               //endpoints of inner integral 
     1904         
     1905        summ = 0.0;                     //initialize intergral 
     1906         
     1907        scale = w[0]; 
     1908        Dnn = w[1];                                     //Nearest neighbor distance A 
     1909        gg = w[2];                                      //Paracrystal distortion factor 
     1910        Rad = w[3];                                     //Sphere radius 
     1911        sld = w[4]; 
     1912        sldSolv = w[5]; 
     1913        background = w[6];  
     1914         
     1915        contrast = sld - sldSolv; 
     1916         
     1917        //Volume fraction calculated from lattice symmetry and sphere radius 
     1918        latticeScale = 2.0*(4.0/3.0)*Pi*(Rad*Rad*Rad)/pow(Dnn/sqrt(3.0/4.0),3); 
     1919         
     1920        for(i=0;i<nordi;i++) { 
     1921                //setup inner integral over the ellipsoidal cross-section 
     1922                summj=0.0; 
     1923                zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;            //the outer dummy is phi 
     1924                for(j=0;j<nordj;j++) { 
     1925                        //20 gauss points for the inner integral 
     1926                        zij = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;               //the inner dummy is theta 
     1927                        yyy = Gauss150Wt[j] * BCC_Integrand(w,x,zi,zij); 
     1928                        summj += yyy; 
     1929                } 
     1930                //now calculate the value of the inner integral 
     1931                answer = (vbj-vaj)/2.0*summj; 
     1932                 
     1933                //now calculate outer integral 
     1934                yyy = Gauss150Wt[i] * answer; 
     1935                summ += yyy; 
     1936        }               //final scaling is done at the end of the function, after the NT_FP64 case 
     1937         
     1938        answer = (vb-va)/2.0*summ; 
     1939        // Multiply by contrast^2 
     1940        answer *= SphereForm_Paracrystal(Rad,contrast,x)*scale*latticeScale; 
     1941        // add in the background 
     1942        answer += background; 
     1943         
     1944        return answer; 
     1945} 
     1946 
     1947// xx is phi (outer) 
     1948// yy is theta (inner) 
     1949double 
     1950BCC_Integrand(double w[], double qq, double xx, double yy) { 
     1951         
     1952        double retVal,temp1,temp3,aa,Da,Dnn,gg,Pi; 
     1953         
     1954        Dnn = w[1]; //Nearest neighbor distance A 
     1955        gg = w[2]; //Paracrystal distortion factor 
     1956        aa = Dnn; 
     1957        Da = gg*aa; 
     1958         
     1959        Pi = 4.0*atan(1.0); 
     1960        temp1 = qq*qq*Da*Da; 
     1961        temp3 = qq*aa;   
     1962         
     1963        retVal = BCCeval(yy,xx,temp1,temp3); 
     1964        retVal /=4.0*Pi; 
     1965         
     1966        return(retVal); 
     1967} 
     1968 
     1969double 
     1970BCCeval(double Theta, double Phi, double temp1, double temp3) { 
     1971 
     1972        double temp6,temp7,temp8,temp9,temp10; 
     1973        double result; 
     1974         
     1975        temp6 = sin(Theta); 
     1976        temp7 = sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)+cos(Theta); 
     1977        temp8 = -1.0*sin(Theta)*cos(Phi)-sin(Theta)*sin(Phi)+cos(Theta); 
     1978        temp9 = -1.0*sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)-cos(Theta); 
     1979        temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9))); 
     1980        result = pow(1.0-(temp10*temp10),3)*temp6/((1.0-2.0*temp10*cos(0.5*temp3*(temp7))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp8))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp9))+(temp10*temp10))); 
     1981         
     1982        return (result); 
     1983} 
     1984 
     1985double 
     1986SphereForm_Paracrystal(double radius, double delrho, double x) {                                         
     1987         
     1988        double bes,f,vol,f2,pi; 
     1989        pi = 4.0*atan(1.0); 
     1990        // 
     1991        //handle q==0 separately 
     1992        if(x==0) { 
     1993                f = 4.0/3.0*pi*radius*radius*radius*delrho*delrho*1.0e8; 
     1994                return(f); 
     1995        } 
     1996         
     1997        bes = 3.0*(sin(x*radius)-x*radius*cos(x*radius))/(x*x*x)/(radius*radius*radius); 
     1998        vol = 4.0*pi/3.0*radius*radius*radius; 
     1999        f = vol*bes*delrho      ;       // [=]  
     2000        // normalize to single particle volume, convert to 1/cm 
     2001        f2 = f * f / vol * 1.0e8;               // [=] 1/cm 
     2002         
     2003        return (f2); 
     2004} 
     2005 
     2006/*      FCC_ParaCrystal  :  calculates the form factor of a Triaxial Ellipsoid at the given x-value p->x 
     2007 
     2008Uses 150 pt Gaussian quadrature for both integrals 
     2009 
     2010*/ 
     2011double 
     2012FCC_ParaCrystal(double w[], double x) 
     2013{ 
     2014        int i,j; 
     2015        double Pi; 
     2016        double scale,Dnn,gg,Rad,contrast,background,latticeScale,sld,sldSolv;           //local variables of coefficient wave 
     2017        int nordi=150;                  //order of integration 
     2018        int nordj=150; 
     2019        double va,vb;           //upper and lower integration limits 
     2020        double summ,zi,yyy,answer;                      //running tally of integration 
     2021        double summj,vaj,vbj,zij;                       //for the inner integration 
     2022         
     2023        Pi = 4.0*atan(1.0); 
     2024        va = 0.0; 
     2025        vb = 2.0*Pi;            //orintational average, outer integral 
     2026        vaj = 0.0; 
     2027        vbj = Pi;               //endpoints of inner integral 
     2028         
     2029        summ = 0.0;                     //initialize intergral 
     2030         
     2031        scale = w[0]; 
     2032        Dnn = w[1];                                     //Nearest neighbor distance A 
     2033        gg = w[2];                                      //Paracrystal distortion factor 
     2034        Rad = w[3];                                     //Sphere radius 
     2035        sld = w[4]; 
     2036        sldSolv = w[5]; 
     2037        background = w[6];  
     2038         
     2039        contrast = sld - sldSolv; 
     2040        //Volume fraction calculated from lattice symmetry and sphere radius 
     2041        latticeScale = 4.0*(4.0/3.0)*Pi*(Rad*Rad*Rad)/pow(Dnn*sqrt(2.0),3); 
     2042         
     2043        for(i=0;i<nordi;i++) { 
     2044                //setup inner integral over the ellipsoidal cross-section 
     2045                summj=0.0; 
     2046                zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;            //the outer dummy is phi 
     2047                for(j=0;j<nordj;j++) { 
     2048                        //20 gauss points for the inner integral 
     2049                        zij = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;               //the inner dummy is theta 
     2050                        yyy = Gauss150Wt[j] * FCC_Integrand(w,x,zi,zij); 
     2051                        summj += yyy; 
     2052                } 
     2053                //now calculate the value of the inner integral 
     2054                answer = (vbj-vaj)/2.0*summj; 
     2055                 
     2056                //now calculate outer integral 
     2057                yyy = Gauss150Wt[i] * answer; 
     2058                summ += yyy; 
     2059        }               //final scaling is done at the end of the function, after the NT_FP64 case 
     2060         
     2061        answer = (vb-va)/2.0*summ; 
     2062        // Multiply by contrast^2 
     2063        answer *= SphereForm_Paracrystal(Rad,contrast,x)*scale*latticeScale; 
     2064        // add in the background 
     2065        answer += background; 
     2066         
     2067        return answer; 
     2068} 
     2069 
     2070 
     2071// xx is phi (outer) 
     2072// yy is theta (inner) 
     2073double 
     2074FCC_Integrand(double w[], double qq, double xx, double yy) { 
     2075         
     2076        double retVal,temp1,temp3,aa,Da,Dnn,gg,Pi; 
     2077         
     2078        Pi = 4.0*atan(1.0); 
     2079        Dnn = w[1]; //Nearest neighbor distance A 
     2080        gg = w[2]; //Paracrystal distortion factor 
     2081        aa = Dnn; 
     2082        Da = gg*aa; 
     2083         
     2084        temp1 = qq*qq*Da*Da; 
     2085        temp3 = qq*aa; 
     2086         
     2087        retVal = FCCeval(yy,xx,temp1,temp3); 
     2088        retVal /=4*Pi; 
     2089         
     2090        return(retVal); 
     2091} 
     2092 
     2093double 
     2094FCCeval(double Theta, double Phi, double temp1, double temp3) { 
     2095 
     2096        double temp6,temp7,temp8,temp9,temp10; 
     2097        double result; 
     2098         
     2099        temp6 = sin(Theta); 
     2100        temp7 = sin(Theta)*sin(Phi)+cos(Theta); 
     2101        temp8 = -1.0*sin(Theta)*cos(Phi)+cos(Theta); 
     2102        temp9 = -1.0*sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi); 
     2103        temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9))); 
     2104        result = pow((1.0-(temp10*temp10)),3)*temp6/((1.0-2.0*temp10*cos(0.5*temp3*(temp7))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp8))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp9))+(temp10*temp10))); 
     2105         
     2106        return (result); 
     2107} 
     2108 
     2109 
     2110/*      SC_ParaCrystal  :  calculates the form factor of a Triaxial Ellipsoid at the given x-value p->x 
     2111 
     2112Uses 150 pt Gaussian quadrature for both integrals 
     2113 
     2114*/ 
     2115double 
     2116SC_ParaCrystal(double w[], double x) 
     2117{ 
     2118        int i,j; 
     2119        double Pi; 
     2120        double scale,Dnn,gg,Rad,contrast,background,latticeScale,sld,sldSolv;           //local variables of coefficient wave 
     2121        int nordi=150;                  //order of integration 
     2122        int nordj=150; 
     2123        double va,vb;           //upper and lower integration limits 
     2124        double summ,zi,yyy,answer;                      //running tally of integration 
     2125        double summj,vaj,vbj,zij;                       //for the inner integration 
     2126         
     2127        Pi = 4.0*atan(1.0); 
     2128        va = 0.0; 
     2129        vb = Pi/2.0;            //orintational average, outer integral 
     2130        vaj = 0.0; 
     2131        vbj = Pi/2.0;           //endpoints of inner integral 
     2132         
     2133        summ = 0.0;                     //initialize intergral 
     2134         
     2135        scale = w[0]; 
     2136        Dnn = w[1];                                     //Nearest neighbor distance A 
     2137        gg = w[2];                                      //Paracrystal distortion factor 
     2138        Rad = w[3];                                     //Sphere radius 
     2139        sld = w[4]; 
     2140        sldSolv = w[5]; 
     2141        background = w[6];  
     2142         
     2143        contrast = sld - sldSolv; 
     2144        //Volume fraction calculated from lattice symmetry and sphere radius 
     2145        latticeScale = (4.0/3.0)*Pi*(Rad*Rad*Rad)/pow(Dnn,3); 
     2146         
     2147        for(i=0;i<nordi;i++) { 
     2148                //setup inner integral over the ellipsoidal cross-section 
     2149                summj=0.0; 
     2150                zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;            //the outer dummy is phi 
     2151                for(j=0;j<nordj;j++) { 
     2152                        //20 gauss points for the inner integral 
     2153                        zij = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;               //the inner dummy is theta 
     2154                        yyy = Gauss150Wt[j] * SC_Integrand(w,x,zi,zij); 
     2155                        summj += yyy; 
     2156                } 
     2157                //now calculate the value of the inner integral 
     2158                answer = (vbj-vaj)/2.0*summj; 
     2159                 
     2160                //now calculate outer integral 
     2161                yyy = Gauss150Wt[i] * answer; 
     2162                summ += yyy; 
     2163        }               //final scaling is done at the end of the function, after the NT_FP64 case 
     2164         
     2165        answer = (vb-va)/2.0*summ; 
     2166        // Multiply by contrast^2 
     2167        answer *= SphereForm_Paracrystal(Rad,contrast,x)*scale*latticeScale; 
     2168        // add in the background 
     2169        answer += background; 
     2170         
     2171        return answer; 
     2172} 
     2173 
     2174// xx is phi (outer) 
     2175// yy is theta (inner) 
     2176double 
     2177SC_Integrand(double w[], double qq, double xx, double yy) { 
     2178         
     2179        double retVal,temp1,temp2,temp3,temp4,temp5,aa,Da,Dnn,gg,Pi; 
     2180         
     2181        Pi = 4.0*atan(1.0); 
     2182        Dnn = w[1]; //Nearest neighbor distance A 
     2183        gg = w[2]; //Paracrystal distortion factor 
     2184        aa = Dnn; 
     2185        Da = gg*aa; 
     2186         
     2187        temp1 = qq*qq*Da*Da; 
     2188        temp2 = pow( 1.0-exp(-1.0*temp1) ,3); 
     2189        temp3 = qq*aa; 
     2190        temp4 = 2.0*exp(-0.5*temp1); 
     2191        temp5 = exp(-1.0*temp1); 
     2192         
     2193         
     2194        retVal = temp2*SCeval(yy,xx,temp3,temp4,temp5); 
     2195        retVal *= 2.0/Pi; 
     2196         
     2197        return(retVal); 
     2198} 
     2199 
     2200double 
     2201SCeval(double Theta, double Phi, double temp3, double temp4, double temp5) { //Function to calculate integrand values for simple cubic structure 
     2202 
     2203        double temp6,temp7,temp8,temp9; //Theta and phi dependent parts of the equation 
     2204        double result; 
     2205         
     2206        temp6 = sin(Theta); 
     2207        temp7 = -1.0*temp3*sin(Theta)*cos(Phi); 
     2208        temp8 = temp3*sin(Theta)*sin(Phi); 
     2209        temp9 = temp3*cos(Theta); 
     2210        result = temp6/((1.0-temp4*cos((temp7))+temp5)*(1.0-temp4*cos((temp8))+temp5)*(1.0-temp4*cos((temp9))+temp5)); 
     2211         
     2212        return (result); 
     2213} 
     2214 
     2215// scattering from a uniform sphere with a Gaussian size distribution 
     2216// 
     2217double 
     2218FuzzySpheres(double dp[], double q) 
     2219{ 
     2220        double pi,x; 
     2221        double scale,rad,pd,sig,rho,rhos,bkg,delrho,sig_surf,f2,bes,vol,f;              //my local names 
     2222        double va,vb,zi,yy,summ,inten; 
     2223        int nord=20,ii; 
     2224         
     2225        pi = 4.0*atan(1.0); 
     2226        x= q; 
     2227         
     2228        scale=dp[0]; 
     2229        rad=dp[1]; 
     2230        pd=dp[2]; 
     2231        sig=pd*rad; 
     2232        sig_surf = dp[3]; 
     2233        rho=dp[4]; 
     2234        rhos=dp[5]; 
     2235        delrho=rho-rhos; 
     2236        bkg=dp[6]; 
     2237         
     2238                         
     2239        va = -4.0*sig + rad; 
     2240        if (va<0) { 
     2241                va=0;           //to avoid numerical error when  va<0 (-ve q-value) 
     2242        } 
     2243        vb = 4.0*sig +rad; 
     2244         
     2245        summ = 0.0;             // initialize integral 
     2246        for(ii=0;ii<nord;ii+=1) { 
     2247                // calculate Gauss points on integration interval (r-value for evaluation) 
     2248                zi = ( Gauss20Z[ii]*(vb-va) + vb + va )/2.0; 
     2249                // calculate sphere scattering 
     2250                // 
     2251                //handle q==0 separately 
     2252                if (x==0.0) { 
     2253                        f2 = 4.0/3.0*pi*zi*zi*zi*delrho*delrho*1.0e8; 
     2254                        f2 *= exp(-0.5*sig_surf*sig_surf*x*x); 
     2255                        f2 *= exp(-0.5*sig_surf*sig_surf*x*x); 
     2256                } else { 
     2257                        bes = 3.0*(sin(x*zi)-x*zi*cos(x*zi))/(x*x*x)/(zi*zi*zi); 
     2258                        vol = 4.0*pi/3.0*zi*zi*zi; 
     2259                        f = vol*bes*delrho;             // [=] A 
     2260                        f *= exp(-0.5*sig_surf*sig_surf*x*x); 
     2261                        // normalize to single particle volume, convert to 1/cm 
     2262                        f2 = f * f / vol * 1.0e8;               // [=] 1/cm 
     2263                } 
     2264         
     2265                yy = Gauss20Wt[ii] *  Gauss_distr(sig,rad,zi) * f2; 
     2266                yy *= 4.0*pi/3.0*zi*zi*zi;              //un-normalize by current sphere volume 
     2267                 
     2268                summ += yy;             //add to the running total of the quadrature 
     2269                 
     2270                 
     2271        } 
     2272        // calculate value of integral to return 
     2273        inten = (vb-va)/2.0*summ; 
     2274         
     2275        //re-normalize by polydisperse sphere volume 
     2276        inten /= (4.0*pi/3.0*rad*rad*rad)*(1.0+3.0*pd*pd); 
     2277         
     2278        inten *= scale; 
     2279        inten += bkg; 
     2280         
     2281    return(inten);      //scale, and add in the background 
     2282} 
     2283 
     2284 
Note: See TracChangeset for help on using the changeset viewer.