source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Alpha/Tinker/FFT_Debye_Spheres.ipf @ 798

#pragma rtGlobals=1             // Use modern global access method.


// with the Debye Sphere method -- as of Jan 2011
// the discretization seems to work now, and is significantly faster than the old double sum
//
// now I need to XOP-it (done)


// routines to calculate I(q) and p(r) from the 3d arrangement of spheres
// using Debye's method
//




//calculate I(q) given a set of spheres as xyz
//
// this is the REALLY slow way
//
Function CalcIQRfromXYZ(xv,yv,zv,qval,rho,rval,grid)
        Wave xv,yv,zv
        Variable qval,rho,rval,grid
        
        Variable num=numpnts(xv),ii,kk,dik
        Variable Iqr,vol,dum,fQR
        
        //NVAR dCall=dCall
        vol=4*Pi/3*rval*rval*rval
        Iqr=0
        fQR=PhiQR(qval,rval)
        //do i=j sum
        for(ii=0;ii<num;ii+=1)
                dum = rho*vol*fQR
                Iqr += dum*dum
        endfor
        //do i != j (double) sum
        for(ii=0;ii<num;ii+=1)
                for(kk=(ii+1);kk<num;kk+=1)
                        dik=AV_Distance(xv[ii],xv[kk],yv[ii],yv[kk],zv[ii],zv[kk]) * grid
                        //dCall += 1
                        dum = rho*vol*fQR
                        Iqr += 2*dum*dum*sinc(dik*qval)
                endfor
        endfor
                
        return(Iqr)
End

//calculate I(q) given a set of spheres as xyz
//
// this is the REALLY slow way
//
Function CalcIQRfromMat(mat,qval,rho,rval,grid)
        Wave mat
        Variable qval,rho,rval,grid
        
        Variable num,ii,kk,dik
        Variable Iqr,vol,dum,fQR,dCall=0
        
        ParseMatrix3D_rho(mat)
        WAVE xv=x3d
        WAVE yv=y3d
        WAVE zv=z3d
        WAVE rho3d=rho3d
        num=numpnts(xv)
        
        vol=4*Pi/3*rval*rval*rval
        Iqr=0
        fQR=PhiQR(qval,rval)
        //do i=j sum
        for(ii=0;ii<num;ii+=1)
                dum = rho3d[ii]*vol*fQR
                Iqr += dum*dum
        endfor
                
        //do i != j (double) sum
        for(ii=0;ii<num;ii+=1)
                for(kk=(ii+1);kk<num;kk+=1)
                        dik=AV_Distance(xv[ii],xv[kk],yv[ii],yv[kk],zv[ii],zv[kk]) * grid
                        //dCall += 1
                        dum = vol*fQR
                        Iqr += 2*dum*dum*rho3d[ii]*rho3d[kk]*sinc(dik*qval)
                endfor
        endfor
        
//      Print "num, dCall = ", num,dCall
        return(Iqr)
End

//calculate I(q) given a set of spheres as xyz
//
// this is the improved way, doing distance binning
// -- substantially speeded up with two XOPs that take care of the double loops
// binning into the histogram is the bottleneck at this point, but I see no 
// easy way to make this multi-processor aware. Supposedly it should be easy
// to do, but I don't see the fast way to dispatch it to threads without writing four
// functions like in the MC simulation code. Maybe data folders?
//
// some of the work is still done in Igor, but it's the really fast stuff, after checking the timing of each step.
//
// can't use the SLD information here, since one of the major assumptions of this simplification
// is that the SLDs of the spheres are all the same
Function CalcIQRfromMat_bin(qW,iqr,xv,yv,zv,rho,rval,grid)
        Wave qW,iqr,xv,yv,zv
        Variable rho,rval,grid
        
        Variable num,ii,kk,dik,iter,t1
        Variable vol,dum,fQR,qval,F2Q
        Variable dmax,binWidth,Qmax,numBins,binIndex,val
        Variable nthreads

        t1 = ticks
        
        num=numpnts(xv)
        
        //find maximum distance
        
//// I don't have enough space for the NxN matrix once there are 10k or so points.
////    Make/O/D/N=(num,num) distMat
////    distMat = 0
//
//      iter = 0
//      for(ii=0;ii<num;ii+=1)
//              for(kk=(ii+1);kk<num;kk+=1)
////                    dist[ii]=AV_Distance(xv[ii],xv[kk],yv[ii],yv[kk],zv[ii],zv[kk]) * grid
//                      distMat[ii][kk]=AV_Distance(xv[ii],xv[kk],yv[ii],yv[kk],zv[ii],zv[kk]) * grid
//              endfor
//      endfor

//      // so this is a quick way to find the maximum distance. entering a value would be easier
//      Make/O/D/N=(num-1) dist
//      for(ii=0;ii<num;ii+=1)
//              for(kk=(ii+1);kk<num;kk+=1)
//                      dist[ii]=AV_Distance(xv[ii],xv[kk],yv[ii],yv[kk],zv[ii],zv[kk]) * grid
//              endfor
//      endfor
//      
//      dmax = WaveMax(dist)
//      Print "dmax = ",dmax
//      Printf "Finding dmax Igor = %g seconds\r",(ticks-t1)/60.15
//      t1 = ticks
        
        ////// Find maxiumum distance
        // write an XOP that takes the xyz and returns the square of the maximum distance, must sqrt then multiply by grid distance
        
        // the  num < 10000 points is an empirical value for where the overhead of threading is OK
        nthreads=ThreadProcessorCount
        if(nthreads == 1 || num < 10000 )               
                dmax = maxDistanceX(xv,yv,zv,0,numpnts(xv))
//              Print "dmax = ",dmax
//              Printf "Finding dmax XOP = %g seconds\r",(ticks-t1)/60.15
//              t1 = ticks
        else
                dmax = maxDistance_Threaded(xv,yv,zv)
//              Print "dmax = ",dmax
        endif
        
        dmax = sqrt(dmax)
        dmax *= grid
//      Print "dmax = ",dmax

//      Printf "Finding dmax Threaded = %g seconds\r",(ticks-t1)/60.15
//      t1 = ticks

        
        // use a bin width Dmax / 10000 as suggested in Otto's book, pg 160
        //
        Variable ndiv=100000
        binWidth = dmax/ndiv
//      Print "binWidth = ",binWidth
        numBins = ndiv
        Make/O/D/N=(numBins) distBins
        
        distBins = 0
        
//      Printf "Initializing distBins = %g seconds\r",(ticks-t1)/60.15
//      t1 = ticks
        
//      for(ii=0;ii<num;ii+=1)
//              for(kk=(ii+1);kk<num;kk+=1)
//                      val = AV_Distance(xv[ii],xv[kk],yv[ii],yv[kk],zv[ii],zv[kk]) * grid
//                      binIndex = trunc(val/binWidth-0.5)
//                      if(binIndex > numBins -1 )
//                              Print "bad index"
//                      else
//                              distBins[binIndex] += 1
//                      endif
//              endfor
//      endfor


///// threading of the bins

        //write an XOP that takes the xyz and bins, and fills them
        //distBins is returned
        nthreads=ThreadProcessorCount

        if(nthreads == 1)
                binDistanceX(xv, yv, zv, distBins, grid, binWidth,0,numpnts(xv))                //the end point is the numpnts, not n-1
        else
                binDistance_Threaded(xv, yv, zv, distBins, grid, binWidth)
        endif

//      Printf "Binning = %g seconds\r",(ticks-t1)/60.15
//      t1 = ticks
        
//      WaveStats/Q distBins
//      Print "binned values, numSpheres = ",V_avg*V_npnts,num
        
        Duplicate/O distBins dist_at_bin
        for(ii=0;ii<numBins;ii+=1)
                dist_at_bin[ii] =  ii*binWidth
        endfor
        
        
        // remove all of the bins with zero distance, then reset the number of Bins
        RemoveZerosXY(distBins, dist_at_bin)
        numBins = numpnts(distBins)
        
//      Printf "Removing zeroes = %g seconds\r",(ticks-t1)/60.15
//      t1 = ticks
        
        vol=4*Pi/3*rval*rval*rval
        
        Iqr=0
        
        Variable nq = numpnts(qW)

//      Make/O/D/N=(nq) Sij_Bin
        
        for(kk=0;kk<nq;kk+=1)           // loop over the q-values
        
                qval = qW[kk]
                fQR=PhiQR(qval,rval)
                F2Q = fQR*fQR*rho*rho*vol*vol
                
                for(ii=0;ii<numBins;ii+=1)              
                        Iqr[kk] +=  distBins[ii]*sinc(qval*dist_at_bin[ii])
                endfor
                
//              Sij_Bin[kk] = num + 2*Iqr[kk]
                
                Iqr[kk] = F2Q*(num + 2*Iqr[kk])
                        
        endfor
        
//      Printf "Loop over q = %g seconds\r",(ticks-t1)/60.15
        
        return(0)
End

// support up to 8 threads at this time
Function binDistance_Threaded(xv, yv, zv, distBins, grid, binWidth)
        Wave xv, yv, zv, distBins
        Variable grid, binWidth
        
        Variable nthreads,mt,left,right
        Variable ii,num

        nthreads=ThreadProcessorCount

        if(nthreads > 8)
                nthreads = 8
        endif
        
        mt = ThreadGroupCreate(nthreads)
        
        ii=0
        left = 0
        num = numpnts(xv)
        
        for(ii=0;ii<nthreads;ii+=1)
                Duplicate/O distBins $("distBins"+num2str(ii))
                
//  this is an even spreading of the points - not appropriate here for the triangle
//                      Print (ii*num/nthreads),((ii+1)*num/nthreads)
//                      ThreadStart mt,i,Cyl_PolyRadius_T(cw,yw,xw,(ii*num/nthreads),((ii+1)*num/nthreads))

// this splits up the triangle into equal area chunks trapezoid left(ii) -> right(ii)
// be sure that for the last iteration of ii, right = 1 *(num)
                right = 1 - sqrt(1- (ii+1)/nthreads)

//                      Print left,right
                
                if(ii==0)
                        Wave distBins0
                        distBins0 = 0
                        //Print (ii*num/nthreads),((ii+1)*num/nthreads)
                        //ThreadStart mt,ii,binDistance_WF(xv, yv, zv, distBins0, grid, binWidth,(ii*num/nthreads),((ii+1)*num/nthreads))
                        ThreadStart mt,ii,binDistance_WF(xv, yv, zv, distBins0, grid, binWidth, left*num, right*num)
                endif
                if(ii==1)
                        Wave distBins1
                        distBins1 = 0
                        ThreadStart mt,ii,binDistance_WF(xv, yv, zv, distBins1, grid, binWidth, left*num, right*num)
                endif
                if(ii==2)
                        Wave distBins2
                        distBins2 = 0
                        ThreadStart mt,ii,binDistance_WF(xv, yv, zv, distBins2, grid, binWidth, left*num, right*num)
                endif
                if(ii==3)
                        Wave distBins3
                        distBins3 = 0
                        ThreadStart mt,ii,binDistance_WF(xv, yv, zv, distBins3, grid, binWidth, left*num, right*num)
                endif
                if(ii==4)
                        Wave distBins4
                        distBins4 = 0
                        ThreadStart mt,ii,binDistance_WF(xv, yv, zv, distBins4, grid, binWidth, left*num, right*num)
                endif
                if(ii==5)
                        Wave distBins5
                        distBins5 = 0
                        ThreadStart mt,ii,binDistance_WF(xv, yv, zv, distBins5, grid, binWidth, left*num, right*num)
                endif
                if(ii==6)
                        Wave distBins6
                        distBins6 = 0
                        ThreadStart mt,ii,binDistance_WF(xv, yv, zv, distBins6, grid, binWidth, left*num, right*num)
                endif
                if(ii==7)
                        Wave distBins7
                        distBins7 = 0
                        ThreadStart mt,ii,binDistance_WF(xv, yv, zv, distBins7, grid, binWidth, left*num, right*num)
                endif
                
                left = right

        endfor
        
        // wait until done
        do
                variable tgs= ThreadGroupWait(mt,100)
        while( tgs != 0 )
        variable dummy= ThreadGroupRelease(mt)
        mt=0
//      Print "done with all threads"

        //      then add them all back together
        if(nthreads == 1)
                distBins = distBins0            // add up each instance
        endif
        if(nthreads == 2)
                distBins = distBins0    + distBins1
        endif
        if(nthreads == 3)
                distBins = distBins0    + distBins1 + distBins2
        endif
        if(nthreads == 4)
                distBins = distBins0    + distBins1 + distBins2 + distBins3
        endif
        if(nthreads == 5)
                distBins = distBins0    + distBins1 + distBins2 + distBins3 + distBins4
        endif
        if(nthreads == 6)
                distBins = distBins0    + distBins1 + distBins2 + distBins3 + distBins4 + distBins5
        endif
        if(nthreads == 7)
                distBins = distBins0    + distBins1 + distBins2 + distBins3 + distBins4 + distBins5 + distBins6
        endif
        if(nthreads == 8)
                distBins = distBins0    + distBins1 + distBins2 + distBins3 + distBins4 + distBins5 + distBins6 + distBins7
        endif
        
        // then clean up
        KillWaves/Z distBins0,distBins1,distBins2,distBins3
        KillWaves/Z distBins4,distBins5,distBins6,distBins7
        
        return(0)
End

// this is just a worker function to get the ThreadStart operation to compile
ThreadSafe Function binDistance_WF(xv, yv, zv, distBins, grid, binWidth,p1,p2)
        Wave xv, yv, zv, distBins
        Variable grid, binWidth,p1,p2

        Variable ret
        ret = binDistanceX(xv, yv, zv, distBins, grid, binWidth,p1,p2)
        
        return(0)
End


// support up to 8 threads at this time
Function maxDistance_Threaded(xv, yv, zv)
        Wave xv, yv, zv
        
        Variable nthreads,mt,left,right
        Variable ii,num

        nthreads=ThreadProcessorCount

        if(nthreads > 8)
                nthreads = 8
        endif
        
        mt = ThreadGroupCreate(nthreads)
        
        ii=0
        left = 0
        num = numpnts(xv)
        
        for(ii=0;ii<nthreads;ii+=1)
                
//  this is an even spreading of the points - not appropriate here for the triangle
//                      Print (ii*num/nthreads),((ii+1)*num/nthreads)
//                      ThreadStart mt,i,Cyl_PolyRadius_T(cw,yw,xw,(ii*num/nthreads),((ii+1)*num/nthreads))

// this splits up the triangle into equal area chunks trapezoid left(ii) -> right(ii)
// be sure that for the last iteration of ii, right = 1 *(num)
                right = 1 - sqrt(1- (ii+1)/nthreads)

//                      Print left,right
                
                ThreadStart mt,ii,maxDistance_WF(xv, yv, zv, left*num, right*num)
                
                left = right

        endfor
        
        // wait until done
        do
                variable tgs= ThreadGroupWait(mt,100)
        while( tgs != 0 )


//      Get the return values, and find the maximum to return
        Variable maxValue = 0
        
        for(ii=0;ii<nthreads;ii+=1)
//                      Print ThreadReturnValue(mt,ii)
                        maxValue = max(maxValue,ThreadReturnValue(mt,ii))
        endfor

        // now release the threads
        variable dummy= ThreadGroupRelease(mt)
        mt=0
//      Print "done with all threads"
        
        return(maxValue)
End

// this is just a worker function to get the ThreadStart operation to compile
ThreadSafe Function maxDistance_WF(xv, yv, zv, p1, p2)
        Wave xv, yv, zv
        Variable p1,p2

        Variable ret
        ret = maxDistanceX(xv,yv,zv,p1,p2)
        
        return(ret)
End

// distance separating two xyz points
// (this is one of the big time-consuming steps)
Function AV_Distance(x1,x2,y1,y2,z1,z2)
        Variable x1,x2,y1,y2,z1,z2
        
        Variable retval
        retval=sqrt( (x1-x2)^2 + (y1-y2)^2 + (z1-z2)^2 )
        return(retval)
End

//bessel function...

Function PhiQR(qval,rval)
        Variable qval,rval
        
        Variable retval,qr
        
        qr=qval*rval
        retval=(sin(qr)-qr*cos(qr))/qr/qr/qr
        
        return(3*retval)
End



//called from the FFT method panel
//
// calls the full double sum that can take the SLD wave into account, if filled
// - as a result, it's slower than the binned calculation
//
Proc DoSpheresCalcFFTPanel(num,qMin,qMax)
        Variable num=100,qmin=0.004,qmax=0.5
        
        Variable t1
        String qStr="qval_full",iStr="ival_full"                //default wave names, always overwritten
        Variable grid
        
        grid=root:FFT_T

        Make/O/D/N=(num) $qStr,$iStr
        $qStr = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))         
        
        Variable estTime,nx
        String str = ""

        nx = NonZeroValues(mat)
        
        estTime = EstimatedTime(nx,num,0)               // 0 =  XOP, 1 = no XOP
        sprintf str, "Estimated time for the calculation is %g seconds. Proceed?",estTime
        DoAlert 1,str
        if(V_Flag==1)           //yes, proceed
                t1=ticks
                fDoCalc($qStr,$iStr,grid,0,1)
//              Printf "Elapsed AltiSpheres time = %g seconds\r\r",(ticks-t1)/60.15
        Endif
End


//called from the FFT method panel
//
// in this method, the distances are binned as by Otto Glatter, and has been partially XOPed
//
Proc DoBinnedSpheresCalcFFTPanel(num,qMin,qMax)
        Variable num=100,qmin=0.004,qmax=0.5
        
        Variable t1
        String qStr="qval_XOP",iStr="ival_XOP"          //default wave names, always overwritten
        Variable grid
        
        grid=root:FFT_T

        Make/O/D/N=(num) $qStr,$iStr
        $qStr = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))         
        
        Variable estTime,nx
        String str = ""

        nx = NonZeroValues(mat)
        
        estTime = EstimatedTime(nx,num,2)               // 0 =  XOP, 1 = no XOP, 2 = binned distances
        sprintf str, "Estimated time for the calculation is %g seconds. Proceed?",estTime
        DoAlert 1,str
        if(V_Flag==1)           //yes, proceed
                t1=ticks
                fDoCalc($qStr,$iStr,grid,2,1)
//              Printf "Elapsed AltiSpheres time = %g seconds\r\r",(ticks-t1)/60.15
        Endif
End

//called from the FFT method panel
//
// in this method, the distances are binned as by Otto Glatter, and has been partially XOPed
//
Proc DoBinnedSLDCalcFFTPanel(num,qMin,qMax)
        Variable num=100,qmin=0.004,qmax=0.5
        
        Variable t1
        String qStr="qval_SLD",iStr="ival_SLD"          //default wave names, always overwritten
        Variable grid
        
        grid=root:FFT_T

        Make/O/D/N=(num) $qStr,$iStr
        $qStr = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))         
        
        Variable estTime,nx
        String str = ""

        nx = NonZeroValues(mat)
        
        estTime = EstimatedTime(nx,num,3)               // 0 =  XOP, 1 = no XOP, 2 = binned distances
        sprintf str, "Estimated time for the calculation is %g seconds. Proceed?",estTime
        DoAlert 1,str
        if(V_Flag==1)           //yes, proceed
                t1=ticks
                fDoCalc($qStr,$iStr,grid,3,1)
//              Printf "Elapsed AltiSpheres time = %g seconds\r\r",(ticks-t1)/60.15
        Endif
End



////////////////////

//waves must exist at this point
// ival_(ext), qval_(ext) are what is passed at this point
//
// switch the type of calculation based on a flag
//
// flag = 0 = full calculation, using double sum
// flag = 2 = binned distances
//                      = 3 = binned distances with SLD
//              = 12 = Ken's generation of XYZ, skip my matrix filling and parsing
//
//
Function fDoCalc(qval,ival,grid,flag,verbose)
        Wave qval,ival
        Variable grid,flag,verbose

        Variable t1,val
        WAVE mat = root:mat
        // convert the matrix to XYZ triplet for calculation
        t1=ticks
        
        if(flag == 12 || flag == 13)                    //skip the parsing
        // when using Ken's method, there is no matrix to parse, or I can put 
        // his xyz (non-integer) onto the mat, then parse it back off (works fine, but not smart)       
                Duplicate/O root:xOutW,x3d
                Duplicate/O root:yOutW,y3d
                Duplicate/O root:zOutW,z3d
                Duplicate/O root:sldW,rho3d
                flag -= 10              //switch to the binned method, either 2 or 3
        else
                ParseMatrix3D_rho(mat)
                if(verbose)
                        Printf "ParseMatrix3D time = %g seconds\r",(ticks-t1)/60.15
                endif
        endif
        
        WAVE x3d=x3d
        WAVE y3d=y3d
        WAVE z3d=z3d
        WAVE rho3d=rho3d
        
        Variable num,nx
        num=numpnts(qval)
        nx=numpnts(x3d)
        
        if(verbose)
                Printf "Estimated  calculation time = %g seconds for %d qvals and %d spheres\r",        EstimatedTime(nx,num,flag),num,nx
        endif
        t1=ticks
        
        NVAR Tscale=root:FFT_T
        NVAR solventSLD = root:FFT_SolventSLD
        
        Variable vol,rval,frac
        rval = grid*0.62
        NVAR delRho = root:FFT_delRho
        vol = nx*(Tscale)^3
        
        frac = nx/DimSize(mat,0)^3
        
        switch(flag)    // numeric switch
                case 0:         // full double sum
                        Duplicate/O rho3D tmpRho3d
                        tmpRho3d = rho3d - solventSLD
                        MultiThread     ival = DebyeSpheresX(qval,x3d,y3d,z3d,tmpRho3d,0.62*grid,grid)
                        ival *= delRho*delRho
                        ival /= vol
                        ival *= 1e8
        
                        ival *= frac
                        
                        KillWaves/Z tmpRho3d
                        break                                   
                case 2:         // binned distances
                                        //////improved version uses distance binning, and is now AAO
                        // don't need to duplicate the wave since I'm passing only a single value for the SLD
                        CalcIQRfromMat_bin(qval,ival,x3d,y3d,z3d,rho3D[0]-solventSLD, rval, grid)
                        ival *= delRho*delRho
                        ival /= vol
                        ival *= 1e8
                        
                        ival *= frac
                        
                        break
                case 3:         // binned distances + SLD
                                        //////improved version uses distance binning, and is now AAO
                        Duplicate/O rho3D tmpRho3d
                        tmpRho3d = rho3d - solventSLD
                        
                        CalcIQRfromMat_bin_SLD(qval,ival,x3d,y3d,z3d,tmpRho3d, rval, grid)
                        ival *= delRho*delRho
                        ival /= vol
                        ival *= 1e8
                        
                        ival *= frac
                        
                        KillWaves/Z tmpRho3d
                        break
                default:                                                        // optional default expression executed
                        Print "No case matched in fDoCalc"                              // when no case matches
        endswitch
        
        if(verbose)
                Printf "Elapsed time for %d qvals and %d spheres = %g seconds\r",num,nx,(ticks-t1)/60.15
        endif
//      Duplicate/O ival_XOP ival_Igor,ival_binned
        
//////
//      Printf "Estimated Igor Function time = %g seconds for %d qvals and %d spheres\r",       EstimatedTime(nx,num,1),num,nx
//      t1=ticks
//      ival_Igor = CalcIQRfromMat(mat,qval,1, rval, grid)
//      
////    vol = nx*4*Pi/3*rval*rval*rval
//      ival_Igor *= delRho*delRho
//      ival_Igor /= vol
//      ival_Igor *= 1e8
//      
//      Printf "Igor function time for %d qvals and %d spheres = %g seconds\r",num,nx,(ticks-t1)/60.15
/////// 
        return(0)
End



// Based on the WM procedure
// RemoveNaNsXY(theXWave, theYWave)
//      Removes all Zeroes in an XY pair if the X wave has the value 0.
//      Returns the number of points removed.
Function RemoveZerosXY(theXWave, theYWave)
        Wave theXWave
        Wave theYWave

        Variable p, numPoints, numNaNs
        Variable xval, yval
        
        numNaNs = 0
        p = 0                                                                                   // the loop index
        numPoints = numpnts(theXWave)                   // number of times to loop

        do
                xval = theXWave[p]
                yval = theYWave[p]
                if (xval == 0)          //is xWave == 0?
                        numNaNs += 1
                else                                                                            // if not an outlier
                        theXWave[p - numNaNs] = xval            // copy to input wave
                        theYWave[p - numNaNs] = yval            // copy to input wave
                endif
                p += 1
        while (p < numPoints)
        
        // Truncate the wave
        DeletePoints numPoints-numNaNs, numNaNs, theXWave, theYWave
        
        return(numNaNs)
End


///////// binned Debye, with multiple SLDs
// as in E. Pantos et.al. J. Molec. Struct. 383 (1996) 303-308.
//
//calculate I(q) given a set of spheres as xyz
//
// this is the improved way, doing distance binning
// - in this method, there can be more than one SLD (different from the solvent)
//
// -(done)- substantially speeded up with two XOPs that take care of the double loops
// binning into the histogram is the bottleneck at this point, but I see no 
// easy way to make this multi-processor aware. Supposedly it should be easy
// to do, but I don't see the fast way to dispatch it to threads without writing four
// functions like in the MC simulation code. Maybe data folders?
//
// some of the work is still done in Igor, but it's the really fast stuff, after checking the timing of each step.
//
//----
// The bookkeeping here is a bit of a mess to keep track of the Gij binning, so here are the details. Note 
// that as ugly and as many loops there are in the function, the time spent in these loops
// is practically nothing (much less than a second) compared to the binning.
//
// SLD_id : a vector with the values of the SLD that are not solvent. The point value of this wave
//                              is used to refer to the "material" with the given SLD, such that SLD_id[index] = SLD_Value
// SLDLookup :  vector that does the reverse correlation to pair up the SLD with its index, such that
//                                              SLDLookup[SLD_Value] = index
// psf_id_vec : a vector of the ij pairs, numbered starting from 11, 12, 13, 22, 23, 33, etc.
// psf_id_mat : a matrix such that psf_id_mat[i][j] gives the INDEX of the (i+1)(j+1) PSF pair. For example,
//                                              psf_id_mat[0][2] corresponds to PSF 13, with and index (in psf_id_vec) of 2.
// Nij          :       a vector of the number of Ni type spheres. Index matches the psf_id_vec. Values are zero
//                                      for i != j, and are the number of spheres for i==j
//
// binMatrix : a matrix of /N=(numBins,NumDiffPSF) that contains the Gij binned pairwise based on the SLD
//                              the columns (numDiffPSF) are indexed to the psf_id_vec
//
// The Gijx and Sijx waves are calculated and saved as a function of r or q, where x is the numerical
// index that corresponds to psf_id_vec. Note that the Sijx are not exactly what you expect for a partial
// structure factor. Some conversion is necessary and the details are still to be worked out, but this is 
// all within Igor, and not in the XOP.
//
Function CalcIQRfromMat_bin_SLD(qW,iqr,xv,yv,zv,rho,rval,grid)
        Wave qW,iqr,xv,yv,zv,rho
        Variable rval,grid
        
        Variable num,ii,jj,kk,dik,iter,t1,np,index,nq
        Variable vol,dum,fQR,qval,F2Q
        Variable dmax,binWidth,Qmax,numBins,binIndex,val
        Variable rhoi,rhok,newNumBins
        Variable nthreads

        NVAR    solventSLD = root:FFT_SolventSLD

//      t1 = ticks
        
        num=numpnts(xv)
        
//      Printf "\r\rParsing = %g seconds\r",(ticks-t1)/60.15
//      t1 = ticks
                
        ////// Find maxiumum distance
        // write an XOP that takes the xyz and returns the square of the maximum distance, must sqrt then multiply by grid distance
        nthreads=ThreadProcessorCount
        if(nthreads == 1)
                dmax = maxDistanceX(xv,yv,zv,0,numpnts(xv))
        else
                dmax = maxDistance_Threaded(xv,yv,zv)
        endif
        
//      dmax = maxDistanceX(xv,yv,zv,0,numpnts(xv))
        dmax = sqrt(dmax)
        dmax *= grid
//      Print "dmax = ",dmax

//      Printf "Finding dmax XOP = %g seconds\r",(ticks-t1)/60.15
//      t1 = ticks

// Now since there are multiple SLDs, we need:
// - a count of the number of different SLDs
// - the number of different PSFs
// - make the histogram matrix, keeping a lookup table of the psf_id

// how many different SLDs are there?
        //loop through the rho wave, disregarding the solvent
        Duplicate/O rho tmpSLD
        Make/O/D/N=0 SLD_id
        SLD_id = 0
        np = 0
        do
                val = WaveMin(tmpSLD)
                if(val != solventSLD && val != 1e10)
                        np += 1
                        InsertPoints numpnts(SLD_id),1,SLD_id
                        SLD_id[np-1] = val
                endif
                tmpSLD = (tmpSLD == val) ? 1e10 : tmpSLD                // set the values to 1e10
        while(val != 1e10)
        KillWaves/Z tmpSLD      

        Variable numDiffSLD = numpnts(SLD_id)
        Variable numDiffPSF, PSFIndex
        numDiffPSF = NumDiffSLD*(numDiffSLD + 1)/2

// make a lookup table to get from rho@xyz to the index
        Variable maxSLD = WaveMax(SLD_id)
        Make/O/D/N=(maxSLD+1) SLDlookup
        SLDLookup = NaN
        ii=0
        do
                SLDLookup[SLD_id[ii]] = ii
                ii+=1
        while(ii<numpnts(SLD_id))

//      Printf "Find number of different SLDs = %g seconds\r",(ticks-t1)/60.15
//      t1 = ticks

        // use a bin width Dmax / 10000 as suggested in Otto's book, pg 160
        //
        numBins=100000
        binWidth = dmax/numBins
        
        Make/O/D/N=(numBins,NumDiffPSF) binMatrix
        binMatrix=0
        
        Make/O/D/N=(numDiffSLD,numDiffSLD) psf_id_mat
        Make/O/D/N=(numDiffPSF) psf_id_vec,Nij
        psf_id_mat = 0
        psf_id_vec = 0
        Nij = 0
        
        index = 0
        for(ii=1;ii<=numDiffSLD;ii+=1)
                for(kk=ii;kk<=numDiffSLD;kk+=1)         
                        //Print ii,kk   
                        psf_id_vec[index] = 10*ii+kk
                        psf_id_mat[ii-1][kk-1] = index
                        psf_id_mat[kk-1][ii-1] = index
                        index += 1      
                endfor
        endfor
        
//      Printf "Set up SLD id matrix, making waves = %g seconds\r",(ticks-t1)/60.15
//      t1 = ticks

        for(ii=0;ii<num;ii+=1)
                rhoi = rho[ii]
                // find the Nij. Not sure if these are what I really need...
                PSFIndex = psf_id_mat[ SLDLookup[rhoi] ][ SLDlookup[rhoi] ]
                Nij[PSFIndex] += 1                      //these are just the Ni counts
        endfor
        
//      Printf "get the Ni counts = %g seconds\r",(ticks-t1)/60.15
//      t1 = ticks
        
// this is the double loop that has been XOPed  
//      for(ii=0;ii<num;ii+=1)
//              for(kk=(ii+1);kk<num;kk+=1)
//                      val = AV_Distance(xv[ii],xv[kk],yv[ii],yv[kk],zv[ii],zv[kk]) * grid
//                      binIndex = trunc(val/binWidth-0.5)
//                      if(binIndex > numBins -1 )
//                              Print "bad index"
//                      else
//                              //figure out which PSF the pair distance belongs to
//                              rhoi = rho[ii]
//                              rhok = rho[kk]
//                              PSFIndex = psf_id_mat[ SLDLookup[rhoi] ][ SLDlookup[rhok] ]
//                              binMatrix[binIndex][PSFIndex] += 1
//                              
//                      endif
//              endfor          
//      endfor

// binMatrix is returned

        nthreads=ThreadProcessorCount
        if(nthreads == 1)
                binSLDDistanceX(xv, yv, zv, rho, binMatrix, SLDLookup, psf_id_mat, grid, binWidth,0,numpnts(xv))
        else
                binSLDDistance_SLD_Threaded(xv, yv, zv, rho, binMatrix, SLDLookup, psf_id_mat, grid, binWidth)
        endif
        
//      Printf "binned the distances and SLDs = %g seconds\r",(ticks-t1)/60.15
//      t1 = ticks

// Now for each Gij
        // -- loop over q, calculating Sij

        nq = numpnts(qW)
        
        for(jj=0;jj<numDiffPSF;jj+=1)
        
                Make/O/D/N=(numBins) dist_at_bin, Gij

                for(ii=0;ii<numBins;ii+=1)
                        dist_at_bin[ii] =  ii*binWidth
                endfor
                                
                Gij = binMatrix[p][jj]
                
                // remove all of the bins with zero distance, then reset the number of Bins
                RemoveZerosXY(Gij, dist_at_bin)
                
                Duplicate/O Iqr, Sij
                Sij = 0
                newNumBins = numpnts(Gij)
                                
                for(kk=0;kk<nq;kk+=1)
                        qval = qW[kk]
                        for(ii=0;ii<newNumBins;ii+=1)           
                                Sij[kk] += Gij[ii]*sinc(qval*dist_at_bin[ii])
                        endfor
                        Sij[kk] *= 2
//                      Sij[kk] += Nij[jj]                      // these may not be the proper Nij
                endfor
                
                // keep copies before looping to the next Gij
                Duplicate/O Sij,$("Sij"+num2str(jj))
                Duplicate/O Gij, $("Gij"+num2str(jj))
                Duplicate/O dist_at_bin,$("dist_at_bin"+num2str(jj)) 
        
        endfor

//      Printf "Calculate Sij's and Gij's = %g seconds\r",(ticks-t1)/60.15
        t1 = ticks
        
        // now add up everything to get the intensity
        // use the psf_id_mat[][] to get the Gij index
        
        vol=4*Pi/3*rval*rval*rval
        Iqr = 0
        
        for(kk=0;kk<nq;kk+=1)           // loop over the q-values
        
                qval = qW[kk]
                fQR = PhiQR(qval,rval)

                for(ii=0;ii<numDiffSLD;ii+=1)
                        index = psf_id_mat[ii][ii]
                        
                        F2Q = fQR*fQR*SLD_id[ii]*SLD_id[ii]*vol*vol
                        
                        Iqr[kk] += Nij[index]*F2Q                       //these are just the "ii" numbers
                
                        for(jj=ii;jj<numDiffSLD;jj+=1)  
                                index = psf_id_mat[ii][jj]
                                
                                WAVE Sij = $("root:Sij"+num2str(index))         //switch the Sij
                                
                                Iqr[kk] += fQR*fQR*SLD_id[ii]*SLD_id[jj]*vol*vol*Sij[kk]                                //factor of 2 removed, not sure why
                                
                        endfor
                endfor
                
        endfor
                
//      Printf "Last loop to calculate I(q) = %g seconds\r",(ticks-t1)/60.15

        return (0)
End


// support up to 8 threads at this time
Function binSLDDistance_SLD_Threaded(xv, yv, zv, rho, binMatrix, SLDLookup, psf_id_mat, grid, binWidth)
        Wave xv, yv, zv, rho, binMatrix, SLDLookup, psf_id_mat
        Variable grid, binWidth
        
        Variable nthreads,mt,left,right
        Variable ii,num

        nthreads=ThreadProcessorCount

        if(nthreads > 8)
                nthreads = 8
        endif
        
        mt = ThreadGroupCreate(nthreads)
        
        ii=0
        left = 0
        num = numpnts(xv)
        
        for(ii=0;ii<nthreads;ii+=1)
                Duplicate/O binMatrix $("binMatrix"+num2str(ii))
                
//  this is an even spreading of the points - not appropriate here for the triangle
//                      Print (ii*num/nthreads),((ii+1)*num/nthreads)
//                      ThreadStart mt,i,Cyl_PolyRadius_T(cw,yw,xw,(ii*num/nthreads),((ii+1)*num/nthreads))

// this splits up the triangle into equal area chunks trapezoid left(ii) -> right(ii)
// be sure that for the last iteration of ii, right = 1 *(num)
                right = 1 - sqrt(1- (ii+1)/nthreads)

//                      Print left,right
                
                if(ii==0)
                        Wave binMatrix0
                        binMatrix0 = 0
                        //Print (ii*num/nthreads),((ii+1)*num/nthreads)
                        //ThreadStart mt,ii,binDistance_WF(xv, yv, zv, binMatrix0, grid, binWidth,(ii*num/nthreads),((ii+1)*num/nthreads))
                        ThreadStart mt,ii,binDistance_SLD_WF(xv, yv, zv, rho, binMatrix0, SLDLookup, psf_id_mat, grid, binWidth, left*num, right*num)
                endif
                if(ii==1)
                        Wave binMatrix1
                        binMatrix1 = 0
                        ThreadStart mt,ii,binDistance_SLD_WF(xv, yv, zv, rho, binMatrix1, SLDLookup, psf_id_mat, grid, binWidth, left*num, right*num)
                endif
                if(ii==2)
                        Wave binMatrix2
                        binMatrix2 = 0
                        ThreadStart mt,ii,binDistance_SLD_WF(xv, yv, zv, rho, binMatrix2, SLDLookup, psf_id_mat, grid, binWidth, left*num, right*num)
                endif
                if(ii==3)
                        Wave binMatrix3
                        binMatrix3 = 0
                        ThreadStart mt,ii,binDistance_SLD_WF(xv, yv, zv, rho, binMatrix3, SLDLookup, psf_id_mat, grid, binWidth, left*num, right*num)
                endif
                if(ii==4)
                        Wave binMatrix4
                        binMatrix4 = 0
                        ThreadStart mt,ii,binDistance_SLD_WF(xv, yv, zv, rho, binMatrix4, SLDLookup, psf_id_mat, grid, binWidth, left*num, right*num)
                endif
                if(ii==5)
                        Wave binMatrix5
                        binMatrix5 = 0
                        ThreadStart mt,ii,binDistance_SLD_WF(xv, yv, zv, rho, binMatrix5, SLDLookup, psf_id_mat, grid, binWidth, left*num, right*num)
                endif
                if(ii==6)
                        Wave binMatrix6
                        binMatrix6 = 0
                        ThreadStart mt,ii,binDistance_SLD_WF(xv, yv, zv, rho, binMatrix6, SLDLookup, psf_id_mat, grid, binWidth, left*num, right*num)
                endif
                if(ii==7)
                        Wave binMatrix7
                        binMatrix7 = 0
                        ThreadStart mt,ii,binDistance_SLD_WF(xv, yv, zv, rho, binMatrix7, SLDLookup, psf_id_mat, grid, binWidth, left*num, right*num)
                endif
                
                left = right

        endfor
        
        // wait until done
        do
                variable tgs= ThreadGroupWait(mt,100)
        while( tgs != 0 )
        variable dummy= ThreadGroupRelease(mt)
        mt=0
//      Print "done with all threads"

        //      then add them all back together
        if(nthreads == 1)
                binMatrix = binMatrix0          // add up each instance
        endif
        if(nthreads == 2)
                binMatrix = binMatrix0  + binMatrix1
        endif
        if(nthreads == 3)
                binMatrix = binMatrix0  + binMatrix1 + binMatrix2
        endif
        if(nthreads == 4)
                binMatrix = binMatrix0  + binMatrix1 + binMatrix2 + binMatrix3
        endif
        if(nthreads == 5)
                binMatrix = binMatrix0  + binMatrix1 + binMatrix2 + binMatrix3 + binMatrix4
        endif
        if(nthreads == 6)
                binMatrix = binMatrix0  + binMatrix1 + binMatrix2 + binMatrix3 + binMatrix4 + binMatrix5
        endif
        if(nthreads == 7)
                binMatrix = binMatrix0  + binMatrix1 + binMatrix2 + binMatrix3 + binMatrix4 + binMatrix5 + binMatrix6
        endif
        if(nthreads == 8)
                binMatrix = binMatrix0  + binMatrix1 + binMatrix2 + binMatrix3 + binMatrix4 + binMatrix5 + binMatrix6 + binMatrix7
        endif
        
        // then clean up
        KillWaves/Z binMatrix0,binMatrix1,binMatrix2,binMatrix3
        KillWaves/Z binMatrix4,binMatrix5,binMatrix6,binMatrix7
        
        return(0)
End

// this is just a worker function to get the ThreadStart operation to compile
ThreadSafe Function binDistance_SLD_WF(xv, yv, zv, rho, binMatrix, SLDLookup, psf_id_mat, grid, binWidth, p1,p2)
        Wave xv, yv, zv, rho, binMatrix, SLDLookup, psf_id_mat
        Variable grid, binWidth, p1,p2

        Variable ret
        ret = binSLDDistanceX(xv, yv, zv, rho, binMatrix, SLDLookup, psf_id_mat, grid, binWidth, p1,p2)

        return(0)
End