Changeset 934


Ignore:
Timestamp:
May 1, 2020 1:06:35 PM (2 years ago)
Author:
ilavsky
Message:

3D models, seems like the 3D solid may actually work now?

Location:
trunk/User Procedures/Irena
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • trunk/User Procedures/Irena/IR1_Main.ipf

    r933 r934  
    203203        SubMenu "3D Models"
    204204                "Mass Fractal Aggregate", IR3A_MassFractalAggregate()
    205                 //"Two Phase Solids", IR3T_TwoPhaseSystem()
    206                 //"Display 3D data", IR3A_Display3DData()
    207                 //"Import POV or PDB", IR3P_ImportPOVPDB()
     205                "Two Phase Solids", IR3T_TwoPhaseSystem()
     206                "Display 3D data", IR3A_Display3DData()
     207                "Import POV or PDB", IR3P_ImportPOVPDB()
    208208        end
    209209        SubMenu "Anisotropy"
     
    19571957                        if(stringmatch(ba.ctrlName,"OpenGCManuscriptWebPage"))
    19581958                                //doi:10.1007/s11661-009-9950-x
    1959                                 BrowseURL "http://www.jomgateway.net/ArticlePage.aspx?DOI=10.1007/s11661-009-9950-x"
     1959                                BrowseURL "http://dx.doi.org/10.1007/s11661-009-9950-x"
    19601960                        endif
    19611961                        if(stringmatch(ba.ctrlName,"OpenMotofitManuscriptWebPage"))
  • trunk/User Procedures/Irena/IR3_3DModels.ipf

    r893 r934  
    764764        //pick the parameters...
    765765        variable voxelSize = PrimarySize/10
    766         variable NumRSteps = 200
    767766        variable IsoValue = 0.5
    768         variable oversample = 1
    769767        variable Qmin = 0.5 * pi/AggSize
    770768        variable Qmax = 6 * pi/PrimarySize
     
    778776                Wave ThreeDVoxelGram = Wave3DwithPrimary                //thsi is voxelgram
    779777        endif
    780         //if the 3DVocelgram is sufficiently small, oversample to get better data...
    781         if(dimsize(ThreeDVoxelGram,0)*dimsize(ThreeDVoxelGram,1)*dimsize(ThreeDVoxelGram,2)< 80^3)
    782                 oversample = 4
    783         elseif(dimsize(ThreeDVoxelGram,0)*dimsize(ThreeDVoxelGram,1)*dimsize(ThreeDVoxelGram,2)< 150^3)
    784                 oversample = 2
    785         else
    786                 oversample = 1
    787         endif
    788 
    789 
    790 abort
     778
    791779        //Calculate pdf intensity
    792         IR3T_CreatePDF(ThreeDVoxelGram,VoxelSize, NumRSteps, IsoValue, oversample, Qmin, Qmax, NumQSteps)
     780        //TODO: check these are set correctly: VoxelSize, NumRSteps
     781        SetScale /P x, 0, VoxelSize, "A" ,ThreeDVoxelGram
     782        SetScale /P y, 0, VoxelSize, "A" ,ThreeDVoxelGram
     783        SetScale /P z, 0, VoxelSize, "A" ,ThreeDVoxelGram
     784        IR3T_CreatePDFIntensity(ThreeDVoxelGram, IsoValue,  Qmin, Qmax, NumQSteps)
    793785        //append to graph...
    794786        Wave PDFQWv
     
    25532545        Button CreateFolder,pos={250,87},size={100,20},proc=IR3P_POVPDBButtonProc,title="Create Folder", help={"Create Folder for data"}
    25542546
    2555         SetVariable CurrentFolderName,value= root:Packages:POVPDBImport:CurrentFolderName, noedit=1,frame=0
     2547        SetVariable CurrentFolderName,value= root:Packages:POVPDBImport:CurrentFolderName
    25562548        SetVariable CurrentFolderName,pos={15,120},size={280,20},title="Current Folder Name",noproc, help={"Current FOlder name to use"}
    25572549
    25582550        TitleBox FakeLine1 title=" ",fixedSize=1,size={330,3},pos={16,145},frame=0,fColor=(0,0,52224), labelBack=(0,0,52224)
    25592551        TitleBox Info2 title="\Zr120Import data",pos={60,160},frame=0,fstyle=1, fixedSize=1,size={300,20},fColor=(0,0,52224)
    2560         Button Import3DData,pos={30,190},size={150,20},proc=IR3P_POVPDBButtonProc,title="Import 3D Data", help={"Import 3D data in this folder"}
    2561         Button Display3DData,pos={30,220},size={150,20},proc=IR3P_POVPDBButtonProc,title="Display 3D Data", help={"Display 3D data in this folder"}
    2562         Button ImportIntQData,pos={30,250},size={150,20},proc=IR3P_POVPDBButtonProc,title="Import Int/Q Data", help={"Import 1D data in this folder"}
    2563         Button DisplayIntQData,pos={30,280},size={150,20},proc=IR3P_POVPDBButtonProc,title="Display Int/Q Data", help={"Display 1D data in this folder"}
    2564         Button CalculateIntQData,pos={30,310},size={150,20},proc=IR3P_POVPDBButtonProc,title="Calculate Int/Q Data", help={"Calculate 1D data and append"}
     2552
     2553
     2554        SetVariable voxelSize,value= root:packages:POVPDBImport:voxelSize, frame=1
     2555        SetVariable voxelSize,pos={15,190},noproc, help={"WHat is voxel Size of the imported 3D structure? "}
     2556        SetVariable voxelSize title="Voxel Size [A]                ",size={200,17},limits={1,100,1}
     2557
     2558
     2559        Button Import3DData,pos={30,220},size={150,20},proc=IR3P_POVPDBButtonProc,title="Import 3D Data", help={"Import 3D data in this folder"}
     2560        Button Display3DData,pos={30,250},size={150,20},proc=IR3P_POVPDBButtonProc,title="Display 3D Data", help={"Display 3D data in this folder"}
     2561        Button ImportIntQData,pos={30,300},size={150,20},proc=IR3P_POVPDBButtonProc,title="Import Int/Q Data", help={"Import 1D data in this folder"}
     2562        Button DisplayIntQData,pos={30,330},size={150,20},proc=IR3P_POVPDBButtonProc,title="Display Int/Q Data", help={"Display 1D data in this folder"}
     2563        Button CalculateIntQData,pos={30,360},size={150,20},proc=IR3P_POVPDBButtonProc,title="Calculate Int/Q Data", help={"Calculate 1D data and append"}
    25652564
    25662565end
     
    26482647        string/g ListOfStrings
    26492648        //here define the lists of variables and strings needed, separate names by ;...
    2650         ListOfVariables="UseForPOV;UseForPDB;"
     2649        ListOfVariables="UseForPOV;UseForPDB;VoxelSize;"
    26512650        ListOfStrings="NewFolderName;CurrentFolderName;"
    26522651        variable i
     
    26592658        endfor 
    26602659               
     2660        NVAR VoxelSize
     2661        if(VoxelSize<1)
     2662                VoxelSize = 1           //default to 1A
     2663        endif
    26612664        setDataFOlder OldDf
    26622665end
     
    27172720        SVAR NewFolderName=root:packages:POVPDBImport:NewFolderName
    27182721        SVAR CurrentFolderName=root:packages:POVPDBImport:CurrentFolderName
    2719         setDataFOlder root:
    2720         NewDataFOlder/O/S $(PossiblyQuoteName(NewFolderName))
    2721         CurrentFolderName = GetDataFolder(1)
     2722        if(Strlen(NewFolderName)>2)
     2723                setDataFOlder root:
     2724                NewDataFOlder/O/S $(PossiblyQuoteName(NewFolderName))
     2725                CurrentFolderName = GetDataFolder(1)
     2726        else
     2727                Abort "No Folder name exists, type in name first"
     2728        endif
    27222729end
    27232730//******************************************************************************************************************************************************
     
    27562763        SetDataFOlder $(CurrentFolderName)
    27572764        Wave ThreeDVoxelGram = POVVoxelWave
    2758         variable voxelSize = 2
    2759         variable NumRSteps=300
     2765        NVAR voxelSize = root:packages:POVPDBImport:voxelSize
     2766        print "Warning = IR3P_Calculate1DDataFile needs to get VoxelSize fixzed, it is fixerd at 2A"
     2767        //variable NumRSteps=300
    27602768        variable IsoValue = 0.5
    2761         variable oversample = 4
    27622769        variable Qmin = 0.001
    2763         variable Qmax = 0.3
     2770        variable Qmax = 0.6
    27642771        variable NumQSteps = 200
    2765         IR3T_CreatePDF(ThreeDVoxelGram,VoxelSize, NumRSteps, IsoValue, oversample, Qmin, Qmax, NumQSteps)
     2772        setScale/P x, 0, VoxelSize, ThreeDVoxelGram
     2773        setScale/P y, 0, VoxelSize, ThreeDVoxelGram
     2774        setScale/P z, 0, VoxelSize, ThreeDVoxelGram
     2775        IR3T_CreatePDFIntensity(ThreeDVoxelGram, IsoValue, Qmin, Qmax, NumQSteps)
    27662776        Wave PDFQWv
    27672777        Wave PDFIntensityWv
    27682778        Wave Qwave
    27692779        Wave IntWave
     2780        //and this does not work for odd number of rows/columns/...
     2781        //IR3T_CalcAutoCorelIntensity(ThreeDVoxelGram, IsoValue, Qmin, Qmax, NumQSteps)
     2782        //these are autocorrelation calculated intensities...
     2783        //Wave AutoCorIntensityWv
     2784        //Wave AutoCorQWv
    27702785        variable InvarModel=areaXY(PDFQWv, PDFIntensityWv )
    27712786        variable InvarData=areaXY(Qwave, IntWave )
    27722787        PDFIntensityWv*=InvarData/InvarModel
     2788        //InvarModel=areaXY(AutoCorQWv, AutoCorIntensityWv )
     2789        //AutoCorIntensityWv*=InvarData/InvarModel
     2790       
    27732791        DOWIndow POV1DGraph
    27742792        if(V_Flag)
     
    27782796                        AppendToGraph/W=POV1DGraph  PDFIntensityWv vs PDFQWv
    27792797                endif
     2798                //CheckDisplayed /W=POV1DGraph AutoCorIntensityWv
     2799                //if(V_flag==0)
     2800                //      AppendToGraph/W=POV1DGraph  AutoCorIntensityWv vs AutoCorQWv
     2801                //endif
    27802802                ModifyGraph lstyle(PDFIntensityWv)=9,lsize(PDFIntensityWv)=3,rgb(PDFIntensityWv)=(1,16019,65535)
    27812803                ModifyGraph mode(PDFIntensityWv)=4,marker(PDFIntensityWv)=19
    27822804                ModifyGraph msize(PDFIntensityWv)=3
     2805                //ModifyGraph lsize(AutoCorIntensityWv)=3,rgb(AutoCorIntensityWv)=(3,52428,1)
     2806                Legend/C/N=text0/A=MC
    27832807        endif
    27842808end
  • trunk/User Procedures/Irena/IR3_3DSupportFunctions.ipf

    r869 r934  
    2323///******************************************************************************************************************************************
    2424///******************************************************************************************************************************************
    25 
    26 Function IR3T_CalcTwoPntsCorFIntensity()
    27         //to be run after running IR3T_Calc3DTwoPntCorrelation
    28         //calculates intensity vs q vercor for results of that function
    29        
    30         setDataFOlder root:Packages:TwoPhaseSolidModel:
    31        
    32         DoWindow TwoPhaseSystemData
    33         if(V_Flag)
    34                 DoWIndow/F TwoPhaseSystemData
    35                 //thsi should exist:
    36                 Wave    TwoPntDebyePhCorrFnct    = root:Packages:TwoPhaseSolidModel:TwoPntDebyePhCorrFnct                       //= Debye phase correlation function (TwoPntCorrelationsWv - TwoPntCorrelationsWv[0]^2, negative values = 0), x-scaling set as above
    37                 Wave    TwoPntDebyePhCorrRadii = root:Packages:TwoPhaseSolidModel:TwoPntDebyePhCorrRadii                        //= Radii wave for TwoPntDebyePhCorrFnct
    38                 Wave/Z PDFQWv = root:Packages:TwoPhaseSolidModel:PDFQWv
    39                 Wave/Z PDFIntensityWv = root:Packages:TwoPhaseSolidModel:PDFIntensityWv
    40                 Wave ExtrapolatedIntensity = root:Packages:TwoPhaseSolidModel:ExtrapolatedIntensity
    41                 Wave ExtrapolatedQvector = root:Packages:TwoPhaseSolidModel:ExtrapolatedQvector
    42                 Wave OriginalQvector = root:Packages:TwoPhaseSolidModel:OriginalQvector
    43                 Wave OriginalIntensity = root:Packages:TwoPhaseSolidModel:OriginalIntensity
    44                 Wave/Z TheoreticalIntensityDACF = root:Packages:TwoPhaseSolidModel:TheoreticalIntensityDACF
    45                 Wave/Z QvecTheorIntensityDACF=root:Packages:TwoPhaseSolidModel:QvecTheorIntensityDACF
    46                 NVAR Qmin = root:Packages:TwoPhaseSolidModel:LowQExtrapolationStart
    47                 NVAR Qmax = root:Packages:TwoPhaseSolidModel:HighQExtrapolationEnd
    48                 //calculate intensity:
    49                 Duplicate/O ExtrapolatedQvector, TwoPntModelIntensity, TwoPntModelQvec
    50                 TwoPntModelIntensity = IR3T_ConvertDACFToInt(TwoPntDebyePhCorrRadii,TwoPntDebyePhCorrFnct,TwoPntModelQvec)     
    51                 variable InvarModel
    52                 variable InvarData
    53                 //need to renormalzie this together...
    54                 InvarModel=areaXY(TwoPntModelQvec, TwoPntModelIntensity, Qmin, TwoPntModelQvec[numpnts(TwoPntModelQvec)-2] )
    55                 InvarData=areaXY(OriginalQvector, OriginalIntensity, Qmin, TwoPntModelQvec[numpnts(TwoPntModelQvec)-2]  )
    56                         TwoPntModelIntensity*=InvarData/InvarModel
    57                         CheckDisplayed /W=TwoPhaseSystemData TwoPntModelIntensity
    58                         if(V_flag==0)
    59                                 AppendToGraph/W=TwoPhaseSystemData  TwoPntModelIntensity vs TwoPntModelQvec
    60                         endif
    61                         ModifyGraph lstyle(TwoPntModelIntensity)=9,lsize(TwoPntModelIntensity)=3,rgb(TwoPntModelIntensity)=(2,39321,1)
    62                         //ModifyGraph mode(TwoPntModelIntensity)=4,marker(TwoPntModelIntensity)=19
    63                         //ModifyGraph msize(TwoPntModelIntensity)=3
    64         endif
    65 
    66 end
    67 
    68 
    69 ///******************************************************************************************************************************************
    70 ///******************************************************************************************************************************************
    71 
    72 Function IR3T_Calc3DTwoPntCorrelation()
    73 
    74 
    75         setDataFOlder root:Packages:TwoPhaseSolidModel:
     25///******************************************************************************************************************************************
     26///******************************************************************************************************************************************
     27
     28//Function IR3T_Calc3DTwoPntCorrelation()
     29//
     30//
     31//      setDataFOlder root:Packages:TwoPhaseSolidModel:
     32//      //this code calculates Two-point autocorelation function on 3D wave.
     33//      //1. Checks wave for sensibility. Needs wave with mostly 0 in it and 1 for minority phase, if needed, it will switch between 0 and 1 so it evaluates minority phase.
     34//      //2. Calculates autocorelation in p, q, and r directions, at most 100x100 vectors in each direction (3x100*100 calculations).
     35//      //3. Copies scaling to output wave, so if the 3DWave has correct x scaling, data have correct x dimension
     36//      //4. creates Debye phase correlation function and its radii wave
     37//      //these are resulting waves:
     38//      //      TwoPntCorrelationsWv                            = two points correlatiton function, TwoPntCorrelationsWv[0] is volume fraction, x-scaling set based on voxel size
     39//      //              TwoPntDebyePhCorrFnct                           = Debye phase correlation function (TwoPntCorrelationsWv - TwoPntCorrelationsWv[0]^2, negative values = 0), x-scaling set as above
     40//      //              TwoPntDebyePhCorrRadii                  = Radii wave for TwoPntDebyePhCorrFnct
     41//      // 3DWave does not have to have same length side, but this has not been tested at all yet.
     42//     
     43//     
     44//      //result is in the sample folder in TwoPntCorrelationsWv which has x scaling set per scaling of input wave.
     45//      Wave/Z My3DWv = root:Packages:TwoPhaseSolidModel:TwoPhaseSolidMatrix
     46//     
     47//      if(!WaveExists(My3DWv))
     48//              abort
     49//      endif
     50//     
     51//      //Check  My3DWv
     52//      Wavestats/Q My3DWv
     53//      variable VOlumeFraction
     54//      //1. Min=0 , max=1
     55//      if(V_min!=0 || V_max!=1)
     56//              Abort "Wave must contain 0 and 1 ONLY, 1 being minority phase, 0 being majority phase"
     57//      endif
     58//      variable MatrixPhase=0
     59//      if(V_avg>0.49)
     60//              MatrixPhase=1
     61//              VOlumeFraction = 1-V_avg
     62//              Print "Two-point corelation function characterized minority phase. Minority phase is expressed by 0 in provided matrix"
     63//      else
     64//              VOlumeFraction = V_avg
     65//              Print "Two-point corelation function characterized minority phase. Minority phase is expressed by 1 in provided matrix"
     66//      endif
     67////    variable startTicks=ticks
     68//      variable pDim, qdim, rdim
     69//      variable pstep, qstep, rstep
     70//      pDim  = DimSize(My3DWv, 0 )
     71//      qDim  = DimSize(My3DWv, 1 )
     72//      rDim  = DimSize(My3DWv, 2 )
     73//      variable pDelta, qDelta, rDelta //these should be voxel sides, if scaling is used for dimansions.
     74//      string pUnits, qUnits, rUnits
     75//      pUnits = WaveUnits(My3DWv, 0 )
     76//      qUnits = WaveUnits(My3DWv, 1 )
     77//      rUnits = WaveUnits(My3DWv, 2 )
     78//      pDelta = DimDelta(My3DWv, 0 )
     79//      qDelta = DimDelta(My3DWv, 1 )
     80//      rDelta = DimDelta(My3DWv, 2 )
     81//      if((pDim*qDim*rDim) < (48*48*48))
     82//              Abort "This 3D object seems too small to evaluate, minimum dimensions are 50^3"
     83//      endif
     84//      if(pDelta!=qDelta || qDelta!=rDelta)
     85//              Abort "This 3D object seems to have different side scaling - voxel sides. You can only analyze object with cubical voxels."
     86//      endif
     87//      //OK, now we may be able to caculate something...
     88//      pstep = ceil(pDim/100)
     89//      qstep = ceil(qDim/100)
     90//      rstep = ceil(rDim/100)
     91//      variable MaxLength=max(pDim, qDim, rDim)
     92//      //this is wave for results.
     93//      Make/O/N=(MaxLength) TwoPntCorrelationsWv
     94//      variable i, j
     95//
     96//      //row (p index)
     97//      For(i=0;i<qDim;i+=qstep)
     98//              For(j=0;j<rDim;j+=rstep)
     99//                      ImageTransform /G=(i)/P=(j) getRow My3DWv
     100//                      Wave W_ExtractedRow
     101//                      if(MatrixPhase)
     102//                              W_ExtractedRow=!W_ExtractedRow[p]
     103//                      endif
     104//                      redimension/S W_ExtractedRow
     105//                      //MatrixOp/Free RowCorrelated = correlate(W_ExtractedRow,W_ExtractedRow,0)
     106//                      //RowCorrelated/=numpnts(RowCorrelated)
     107//                      Correlate /AUTO W_ExtractedRow, W_ExtractedRow
     108//                      Wave/Z RowCorrelations
     109//                      if(!WaveExists(RowCorrelations))
     110//                              Duplicate/O W_ExtractedRow, RowCorrelations
     111//                      else
     112//                              RowCorrelations+=W_ExtractedRow
     113//                      endif
     114//              endfor
     115//      endfor
     116//      //RowCorrelations/=(round(qDim/qstep)*round(rDim/rstep))
     117//      //circular correlation causes this to be mirrored around center. All we get is half of the distance across this way
     118//      //reverse RowCorrelations
     119//      //redimension/N=(numpnts(RowCorrelations)/2)     RowCorrelations
     120//      DeletePoints 0, (numpnts(RowCorrelations)/2),  RowCorrelations
     121//     
     122//      //column (q index)
     123//      For(i=0;i<pDim;i+=pstep)
     124//              For(j=0;j<rDim;j+=rstep)
     125//                      ImageTransform /G=(i) /P=(j) getCol My3DWv
     126//                      Wave W_ExtractedCol
     127//                      if(MatrixPhase)
     128//                              W_ExtractedCol=!W_ExtractedCol[p]
     129//                      endif
     130//                      redimension/S W_ExtractedCol
     131//                      //MatrixOp/Free ColCorrelated = correlate(W_ExtractedCol,W_ExtractedCol,0)
     132//                      //ColCorrelated/=numpnts(ColCorrelated)
     133//                      Correlate /AUTO W_ExtractedCol, W_ExtractedCol
     134//                      //ColumnCorrelations+=W_ExtractedCol
     135//                      Wave/Z ColumnCorrelations
     136//                      if(!WaveExists(ColumnCorrelations))
     137//                              Duplicate/O W_ExtractedCol, ColumnCorrelations
     138//                      else
     139//                              ColumnCorrelations+=W_ExtractedCol
     140//                      endif
     141//              endfor
     142//      endfor
     143//      //ColumnCorrelations/=(round(pDim/pstep)*round(rDim/rstep))
     144//      //circular correlation causes this to be mirrored around center. ALl we get is half of the distance across this way
     145//      //reverse ColumnCorrelations
     146//      //redimension/N=(numpnts(ColumnCorrelations)/2)  ColumnCorrelations
     147//      DeletePoints 0, (numpnts(ColumnCorrelations)/2),  ColumnCorrelations
     148//
     149//      //beam (r index)       
     150//      For(i=0;i<pDim;i+=pstep)
     151//              For(j=0;j<qDim;j+=qstep)
     152//                      ImageTransform /BEAM={(i),(j)} getBeam My3DWv
     153//                      Wave W_Beam
     154//                      if(MatrixPhase)
     155//                              W_Beam=!W_Beam[p]
     156//                      endif
     157//                      redimension/S W_Beam
     158//                      //MatrixOp/Free BeamCorrelated = correlate(W_Beam,W_Beam,4)
     159//                      Correlate /AUTO W_Beam, W_Beam
     160//                      //BeamCorrelated/=numpnts(BeamCorrelated)
     161//                      //BeamCorrelations+=BeamCorrelated
     162//                      Wave/Z BeamCorrelations
     163//                      if(!WaveExists(BeamCorrelations))
     164//                              Duplicate/O W_Beam, BeamCorrelations
     165//                      else
     166//                              BeamCorrelations+=W_Beam
     167//                      endif
     168//                      //BeamCorrelations+=W_Beam
     169//              endfor
     170//      endfor
     171//      //BeamCorrelations/=(round(pDim/pstep)*round(qDim/qstep))
     172//      //circular correlation causes this to be mirrored around center. All we get is half of the distance across this way
     173//      //reverse BeamCorrelations
     174//      //redimension/N=(numpnts(BeamCorrelations)/2) BeamCorrelations
     175//      //average, this handles wave of different lenghts.     
     176//      //Correlate/AUTO needs rotation around end point
     177//      DeletePoints 0, (numpnts(BeamCorrelations)/2),  BeamCorrelations
     178//
     179//      IR3T_Average3Waves(RowCorrelations,ColumnCorrelations,BeamCorrelations ,TwoPntCorrelationsWv)
     180//      KillWaves/Z RowCorrelations,ColumnCorrelations,BeamCorrelations
     181////    Duplicate/O BeamCorrelations, TwoPntCorrelationsWv
     182//      //redimension/N=(numpnts(TwoPntCorrelationsWv)/2) TwoPntCorrelationsWv
     183//      //normalize to volume fratcion as epxected
     184//      Wavestats/Q TwoPntCorrelationsWv
     185//      TwoPntCorrelationsWv = VOlumeFraction* TwoPntCorrelationsWv[p]/V_max
     186//     
     187//      Duplicate/O TwoPntCorrelationsWv, TwoPntDebyePhCorrFnct, TwoPntDebyePhCorrRadii
     188//      //??? TwoPntDebyePhCorrFnct =  TwoPntCorrelationsWv - VolumeFraction^2                          //this converts this from TwoPointCor FUnction to Debye Phase Correlation function
     189//      TwoPntDebyePhCorrFnct =  TwoPntCorrelationsWv + VolumeFraction^2                                //this converts this from TwoPointCor FUnction to Debye Phase Correlation function
     190//      //is the above negative or positive? I think it is positive, since the random chanse should be volumefraction^2
     191//     
     192//      TwoPntDebyePhCorrFnct = TwoPntDebyePhCorrFnct[p]>0 ? TwoPntDebyePhCorrFnct : 0                  //this needs to be non-negative... - note it should be bit more complicated, but this is OK for now...
     193//
     194//      SetScale/P x 0,(pDelta),pUnits, TwoPntCorrelationsWv, TwoPntDebyePhCorrFnct                             //assign x-scaling in A from 3D model.
     195//      //Seems like this is distance distribution, not radius distribution...
     196//      //that is why we need the 1/2 there???
     197//
     198//      TwoPntDebyePhCorrRadii = pnt2x(TwoPntDebyePhCorrFnct, p )
     199//     
     200//     
     201//      //Display/K=1 TwoPntCorrelationsWv
     202//      //print (ticks-startTicks)/60
     203//end
     204///******************************************************************************************************************************************
     205///******************************************************************************************************************************************
     206///******************************************************************************************************************************************
     207
     208// I M P O R T A N T - verified for spheres on 5-1-2020, seems to work for:
     209// 3D voxelgram, 0 is empty space, 1 is particle
     210// checked on sphere model with voxel sizes 1, 2, and 4 A and it works.
     211
     212Function IR3T_CalcAutoCorelIntensity(My3DWv, IsoValue, Qmin, Qmax, NumQSteps)
     213        Wave My3DWv
     214        variable IsoValue, Qmin, Qmax, NumQSteps
     215        //same as IR3T_CreatePDFIntensity
     216        //              old:   does the same thing as Function IR3T_Calc3DTwoPntCorrelation()
     217        //                                      but uses 3DAutcorrelation
     218       
     219        print "IsoValue value in IR3T_CalcAutoCorelIntensity is for now not used."
     220        setDataFOlder GetWavesDataFolder(My3DWv, 1 )
     221       
    76222        //this code calculates Two-point autocorelation function on 3D wave.
    77223        //1. Checks wave for sensibility. Needs wave with mostly 0 in it and 1 for minority phase, if needed, it will switch between 0 and 1 so it evaluates minority phase.
     
    84230        //              TwoPntDebyePhCorrRadii                  = Radii wave for TwoPntDebyePhCorrFnct
    85231        // 3DWave does not have to have same length side, but this has not been tested at all yet.
    86        
     232        //next it calculates Intensity vs Q vector same as IR3T_CreatePDFIntensity
    87233       
    88234        //result is in the sample folder in TwoPntCorrelationsWv which has x scaling set per scaling of input wave.
    89         Wave/Z My3DWv = root:Packages:TwoPhaseSolidModel:TwoPhaseSolidMatrix
    90        
    91         if(!WaveExists(My3DWv))
    92                 abort
    93         endif
    94        
     235        print "Check input wave and get dimension values. this can take a logng time. 512^3 wave takes ~30 seconds on fast computer"
     236        //variable starttik=ticks
    95237        //Check  My3DWv
    96         Wavestats/Q My3DWv
    97         variable VOlumeFraction
    98         //1. Min=0 , max=1
     238        variable VolumeFraction
     239        Wavestats/Q/M=1 My3DWv
     240        //Check that Min=0 , max=1
    99241        if(V_min!=0 || V_max!=1)
    100242                Abort "Wave must contain 0 and 1 ONLY, 1 being minority phase, 0 being majority phase"
     
    103245        if(V_avg>0.49)
    104246                MatrixPhase=1
    105                 VOlumeFraction = 1-V_avg
     247                VolumeFraction = 1-V_avg
    106248                Print "Two-point corelation function characterized minority phase. Minority phase is expressed by 0 in provided matrix"
    107249        else
    108                 VOlumeFraction = V_avg
     250                VolumeFraction = V_avg
    109251                Print "Two-point corelation function characterized minority phase. Minority phase is expressed by 1 in provided matrix"
    110252        endif
    111 //      variable startTicks=ticks
    112253        variable pDim, qdim, rdim
    113254        variable pstep, qstep, rstep
     
    115256        qDim  = DimSize(My3DWv, 1 )
    116257        rDim  = DimSize(My3DWv, 2 )
    117         variable pDelta, qDelta, rDelta //these should be voxel sides, if scaling is used for dimansions.
     258        variable pDelta, qDelta, rDelta                                         //these should be voxel sides, if scaling is used for dimansions.
    118259        string pUnits, qUnits, rUnits
    119260        pUnits = WaveUnits(My3DWv, 0 )
     
    129270                Abort "This 3D object seems to have different side scaling - voxel sides. You can only analyze object with cubical voxels."
    130271        endif
    131         //OK, now we may be able to conaculate something...
    132         pstep = ceil(pDim/100)
    133         qstep = ceil(qDim/100)
    134         rstep = ceil(rDim/100)
    135         variable MaxLength=max(pDim, qDim, rDim)
    136         //this is wave for results.
    137         Make/O/N=(MaxLength) TwoPntCorrelationsWv
    138         variable i, j
    139 
    140         //row (p index)
    141         For(i=0;i<qDim;i+=qstep)
    142                 For(j=0;j<rDim;j+=rstep)
    143                         ImageTransform /G=(i)/P=(j) getRow My3DWv
    144                         Wave W_ExtractedRow
    145                         if(MatrixPhase)
    146                                 W_ExtractedRow=!W_ExtractedRow[p]
    147                         endif
    148                         redimension/S W_ExtractedRow
    149                         //MatrixOp/Free RowCorrelated = correlate(W_ExtractedRow,W_ExtractedRow,0)
    150                         //RowCorrelated/=numpnts(RowCorrelated)
    151                         Correlate /AUTO W_ExtractedRow, W_ExtractedRow
    152                         Wave/Z RowCorrelations
    153                         if(!WaveExists(RowCorrelations))
    154                                 Duplicate/O W_ExtractedRow, RowCorrelations
    155                         else
    156                                 RowCorrelations+=W_ExtractedRow
    157                         endif
    158                 endfor
    159         endfor
    160         //RowCorrelations/=(round(qDim/qstep)*round(rDim/rstep))
    161         //circular correlation causes this to be mirrored around center. All we get is half of the distance across this way
    162         //reverse RowCorrelations
    163         //redimension/N=(numpnts(RowCorrelations)/2)     RowCorrelations
    164         DeletePoints 0, (numpnts(RowCorrelations)/2),  RowCorrelations
    165        
    166         //column (q index)
    167         For(i=0;i<pDim;i+=pstep)
    168                 For(j=0;j<rDim;j+=rstep)
    169                         ImageTransform /G=(i) /P=(j) getCol My3DWv
    170                         Wave W_ExtractedCol
    171                         if(MatrixPhase)
    172                                 W_ExtractedCol=!W_ExtractedCol[p]
    173                         endif
    174                         redimension/S W_ExtractedCol
    175                         //MatrixOp/Free ColCorrelated = correlate(W_ExtractedCol,W_ExtractedCol,0)
    176                         //ColCorrelated/=numpnts(ColCorrelated)
    177                         Correlate /AUTO W_ExtractedCol, W_ExtractedCol
    178                         //ColumnCorrelations+=W_ExtractedCol
    179                         Wave/Z ColumnCorrelations
    180                         if(!WaveExists(ColumnCorrelations))
    181                                 Duplicate/O W_ExtractedCol, ColumnCorrelations
    182                         else
    183                                 ColumnCorrelations+=W_ExtractedCol
    184                         endif
    185                 endfor
    186         endfor
    187         //ColumnCorrelations/=(round(pDim/pstep)*round(rDim/rstep))
    188         //circular correlation causes this to be mirrored around center. ALl we get is half of the distance across this way
    189         //reverse ColumnCorrelations
    190         //redimension/N=(numpnts(ColumnCorrelations)/2)  ColumnCorrelations
    191         DeletePoints 0, (numpnts(ColumnCorrelations)/2),  ColumnCorrelations
    192 
    193         //beam (r index)       
    194         For(i=0;i<pDim;i+=pstep)
    195                 For(j=0;j<qDim;j+=qstep)
    196                         ImageTransform /BEAM={(i),(j)} getBeam My3DWv
    197                         Wave W_Beam
    198                         if(MatrixPhase)
    199                                 W_Beam=!W_Beam[p]
    200                         endif
    201                         redimension/S W_Beam
    202                         //MatrixOp/Free BeamCorrelated = correlate(W_Beam,W_Beam,4)
    203                         Correlate /AUTO W_Beam, W_Beam
    204                         //BeamCorrelated/=numpnts(BeamCorrelated)
    205                         //BeamCorrelations+=BeamCorrelated
    206                         Wave/Z BeamCorrelations
    207                         if(!WaveExists(BeamCorrelations))
    208                                 Duplicate/O W_Beam, BeamCorrelations
    209                         else
    210                                 BeamCorrelations+=W_Beam
    211                         endif
    212                         //BeamCorrelations+=W_Beam
    213                 endfor
    214         endfor
    215         //BeamCorrelations/=(round(pDim/pstep)*round(qDim/qstep))
    216         //circular correlation causes this to be mirrored around center. All we get is half of the distance across this way
    217         //reverse BeamCorrelations
    218         //redimension/N=(numpnts(BeamCorrelations)/2) BeamCorrelations
    219         //average, this handles wave of different lenghts.     
    220         //Correlate/AUTO needs rotation around end point
    221         DeletePoints 0, (numpnts(BeamCorrelations)/2),  BeamCorrelations
    222 
    223         IR3T_Average3Waves(RowCorrelations,ColumnCorrelations,BeamCorrelations ,TwoPntCorrelationsWv)
    224         KillWaves/Z RowCorrelations,ColumnCorrelations,BeamCorrelations
    225 //      Duplicate/O BeamCorrelations, TwoPntCorrelationsWv
    226         //redimension/N=(numpnts(TwoPntCorrelationsWv)/2) TwoPntCorrelationsWv
    227         //normalize to volume fratcion as epxected
    228         Wavestats/Q TwoPntCorrelationsWv
    229         TwoPntCorrelationsWv = VOlumeFraction* TwoPntCorrelationsWv[p]/V_max
    230        
    231         Duplicate/O TwoPntCorrelationsWv, TwoPntDebyePhCorrFnct, TwoPntDebyePhCorrRadii
    232         TwoPntDebyePhCorrFnct =  TwoPntCorrelationsWv - VolumeFraction^2                                //thsi converts this from TwoPointCor FUnction to Debye Phase Correlation function
    233        
    234         TwoPntDebyePhCorrFnct = TwoPntDebyePhCorrFnct[p]>0 ? TwoPntDebyePhCorrFnct : 0                  //this needs to be non-negative... - note it should be bit more complicated, but this is OK for now...
    235 
    236         SetScale/P x 0,(pDelta),pUnits, TwoPntCorrelationsWv, TwoPntDebyePhCorrFnct                             //assign x-scaling in A from 3D model.
    237         //Seems like this is distance distribution, not radius distribution...
    238         //that is why we need the 1/2 there???
    239 
     272        Print "Preparation phase done, starting autocorrelation. This may be slow. "
     273        //calculate autocorelation
     274        IR3T_Autocorelate3D(My3DWv)
     275        wave AutoCorMatrix              //this is resulting autcorrelation wave
     276        Print "Finished autocorrelation. Calculate Two Point Corelation function now."
     277        // Extract radial profile...
     278        IR3T_CalcRadialAveProfile3D(AutoCorMatrix)     
     279        KillWaves/Z AutoCorMatrix
     280        Wave RadialWaveProfile                                                                                                  //radial profile with the x voxel scaling...
     281        SetScale /P x, 0, pDelta, "A", RadialWaveProfile                                        //here we set x scaling to match input scaling of the 3D wave voxel size.
     282        Print "Finished Two Point Correlation function calculation. Finish by scaling."
     283        Duplicate/free  RadialWaveProfile, TwoPntCorrelationsWv, TwoPntDebyePhCorrFnct, TwoPntDebyePhCorrRadii
     284        KillWaves/Z RadialWaveProfile
     285        //RadialWaveProfile[0] = 1 and at infinity is = VolumeFraction.
     286        //scaling to match volume fraction needed,
     287        TwoPntCorrelationsWv = VOlumeFraction * TwoPntCorrelationsWv - VOlumeFraction^2
     288
     289        FindLevel/EDGE=2/Q/P TwoPntCorrelationsWv, 0
     290        if(V_Flag==0)
     291                //Duplicate/Free TwoPntCorrelationsWv, PDFWave, PrWave
     292                //PrWave = RadialWaveProfile * x
     293                //PDFWave = PrWave * x^2
     294                //variable CalculatedRg=area(PDFWave,0,V_LevelX)/area(PrWave,0,V_LevelX)
     295                //print "Calculated Rg [A] = "+num2str(CalculatedRg)
     296                TwoPntCorrelationsWv[V_LevelX+1, ] = 0
     297        endif
     298
     299        //the commented part below screwed up the results tremendously. Which is bit weird. In any case, TwoPntDebyePhCorrFnct = TwoPntCorrelationsWv works just fine.
     300        TwoPntDebyePhCorrFnct = TwoPntCorrelationsWv                            // (VolumeFraction+VolumeFraction^2)*TwoPntCorrelationsWv + VolumeFraction^2                            //this converts this from TwoPointCor Function to Debye Phase Correlation function
     301                                                                                                                                                        //or is it: TwoPntDebyePhCorrFnct =  TwoPntCorrelationsWv - VolumeFraction^2    ???
    240302        TwoPntDebyePhCorrRadii = pnt2x(TwoPntDebyePhCorrFnct, p )
    241        
    242        
    243         //Display/K=1 TwoPntCorrelationsWv
    244         //print (ticks-startTicks)/60
    245 end
    246 ///******************************************************************************************************************************************
    247 ///******************************************************************************************************************************************
    248 ///******************************************************************************************************************************************
    249 Function IR3T_Calc3D_3DTwoPntCorrelation()
    250         //does the asme thing as Function IR3T_Calc3DTwoPntCorrelation()
    251         //but uses 3DAutcorrelation
    252 
    253 
    254         setDataFOlder root:Packages:TwoPhaseSolidModel:
    255         //this code calculates Two-point autocorelation function on 3D wave.
    256         //1. Checks wave for sensibility. Needs wave with mostly 0 in it and 1 for minority phase, if needed, it will switch between 0 and 1 so it evaluates minority phase.
    257         //2. Calculates autocorelation in p, q, and r directions, at most 100x100 vectors in each direction (3x100*100 calculations).
    258         //3. Copies scaling to output wave, so if the 3DWave has correct x scaling, data have correct x dimension
    259         //4. creates Debye phase correlation function and its radii wave
    260         //these are resulting waves:
    261         //      TwoPntCorrelationsWv                            = two points correlatiton function, TwoPntCorrelationsWv[0] is volume fraction, x-scaling set based on voxel size
    262         //              TwoPntDebyePhCorrFnct                           = Debye phase correlation function (TwoPntCorrelationsWv - TwoPntCorrelationsWv[0]^2, negative values = 0), x-scaling set as above
    263         //              TwoPntDebyePhCorrRadii                  = Radii wave for TwoPntDebyePhCorrFnct
    264         // 3DWave does not have to have same length side, but this has not been tested at all yet.
    265        
    266         //result is in the sample folder in TwoPntCorrelationsWv which has x scaling set per scaling of input wave.
    267         Wave/Z My3DWv = root:Packages:TwoPhaseSolidModel:TwoPhaseSolidMatrix
    268        
    269         if(!WaveExists(My3DWv))
    270                 abort
    271         endif
    272        
    273         //Check  My3DWv
    274         Wavestats/Q My3DWv
    275         variable VOlumeFraction
    276         //1. Min=0 , max=1
    277         if(V_min!=0 || V_max!=1)
    278                 Abort "Wave must contain 0 and 1 ONLY, 1 being minority phase, 0 being majority phase"
    279         endif
    280         variable MatrixPhase=0
    281         if(V_avg>0.49)
    282                 MatrixPhase=1
    283                 VOlumeFraction = 1-V_avg
    284                 Print "Two-point corelation function characterized minority phase. Minority phase is expressed by 0 in provided matrix"
    285         else
    286                 VOlumeFraction = V_avg
    287                 Print "Two-point corelation function characterized minority phase. Minority phase is expressed by 1 in provided matrix"
    288         endif
    289 //      variable startTicks=ticks
    290         variable pDim, qdim, rdim
    291         variable pstep, qstep, rstep
    292         pDim  = DimSize(My3DWv, 0 )
    293         qDim  = DimSize(My3DWv, 1 )
    294         rDim  = DimSize(My3DWv, 2 )
    295         variable pDelta, qDelta, rDelta //these should be voxel sides, if scaling is used for dimansions.
    296         string pUnits, qUnits, rUnits
    297         pUnits = WaveUnits(My3DWv, 0 )
    298         qUnits = WaveUnits(My3DWv, 1 )
    299         rUnits = WaveUnits(My3DWv, 2 )
    300         pDelta = DimDelta(My3DWv, 0 )
    301         qDelta = DimDelta(My3DWv, 1 )
    302         rDelta = DimDelta(My3DWv, 2 )
    303         if((pDim*qDim*rDim) < (48*48*48))
    304                 Abort "This 3D object seems too small to evaluate, minimum dimensions are 50^3"
    305         endif
    306         if(pDelta!=qDelta || qDelta!=rDelta)
    307                 Abort "This 3D object seems to have different side scaling - voxel sides. You can only analyze object with cubical voxels."
    308         endif
    309         //OK, now we may be able to conaculate something...
    310 
    311         IR3T_Autocorelate3D(My3DWv)
    312        
    313         wave AutoCorMatrix              //this is resulting autcorrelation wave
    314        
    315         //we need better method, but for now we can simply extract beam along z direction and get the G(r)
    316        
    317         ImageTransform /BEAM={(DimSize(AutoCorMatrix, 0)/2 -1),(DimSize(AutoCorMatrix, 1)/2 -1)} getBeam AutoCorMatrix
    318         Wave W_Beam
    319         //this has proper scaling from -Size to +Size along that direction.
    320        
    321        
    322 end
     303        //print (ticks-starttik)/60
     304       
     305        make/O/N=(NumQSteps)/D  AutoCorIntensityWv, AutoCorQWv
     306        AutoCorQWv =    Qmin + p*(Qmax-Qmin)/(NumQSteps-1) 
     307        IN2G_ConvertTologspacing(AutoCorQWv,0)                                                                  //creates log-q spacing in the PDFQWv
     308        //multithread AutoCorIntensityWv =  IR3T_CalcIntensityPDF(PDFQWv[p],PDFWave,RadiiWave)                                                                          //this is PDF from IR3T_CreatePDFIntensity
     309        multithread AutoCorIntensityWv = IR3T_ConvertDACFToInt(TwoPntDebyePhCorrRadii,TwoPntDebyePhCorrFnct,AutoCorQWv)         //and this is equivalent using Autocoreelation function
     310
     311end
     312//Function IR3T_CalcTwoPntsCorFIntensity()
     313//      //to be run after running IR3T_Calc3D_3DTwoPntCorrelation
     314//      //calculates intensity vs q vercor for results of that function
     315//     
     316//      setDataFOlder root:Packages:TwoPhaseSolidModel:
     317//     
     318//      DoWindow TwoPhaseSystemData
     319//      if(V_Flag)
     320//              DoWIndow/F TwoPhaseSystemData
     321//              //thsi should exist:
     322//              Wave    TwoPntDebyePhCorrFnct    = root:Packages:TwoPhaseSolidModel:TwoPntDebyePhCorrFnct                       //= Debye phase correlation function (TwoPntCorrelationsWv - TwoPntCorrelationsWv[0]^2, negative values = 0), x-scaling set as above
     323//              Wave    TwoPntDebyePhCorrRadii = root:Packages:TwoPhaseSolidModel:TwoPntDebyePhCorrRadii                        //= Radii wave for TwoPntDebyePhCorrFnct
     324//              Wave/Z PDFQWv = root:Packages:TwoPhaseSolidModel:PDFQWv
     325//              Wave/Z PDFIntensityWv = root:Packages:TwoPhaseSolidModel:PDFIntensityWv
     326//              Wave ExtrapolatedIntensity = root:Packages:TwoPhaseSolidModel:ExtrapolatedIntensity
     327//              Wave ExtrapolatedQvector = root:Packages:TwoPhaseSolidModel:ExtrapolatedQvector
     328//              Wave OriginalQvector = root:Packages:TwoPhaseSolidModel:OriginalQvector
     329//              Wave OriginalIntensity = root:Packages:TwoPhaseSolidModel:OriginalIntensity
     330//              Wave/Z TheoreticalIntensityDACF = root:Packages:TwoPhaseSolidModel:TheoreticalIntensityDACF
     331//              Wave/Z QvecTheorIntensityDACF=root:Packages:TwoPhaseSolidModel:QvecTheorIntensityDACF
     332//              NVAR Qmin = root:Packages:TwoPhaseSolidModel:LowQExtrapolationStart
     333//              NVAR Qmax = root:Packages:TwoPhaseSolidModel:HighQExtrapolationEnd
     334//              //calculate intensity:
     335//              Duplicate/O ExtrapolatedQvector, TwoPntModelIntensity, TwoPntModelQvec
     336//              TwoPntModelIntensity = IR3T_ConvertDACFToInt(TwoPntDebyePhCorrRadii,TwoPntDebyePhCorrFnct,TwoPntModelQvec)     
     337//              //to make Guinier here match, we need to do this... WHy?
     338//      //      TwoPntModelQvec*=pi
     339//              //end of why???
     340//              variable InvarModel
     341//              variable InvarData
     342//              //need to renormalize this together...
     343//              InvarModel=areaXY(TwoPntModelQvec, TwoPntModelIntensity, Qmin, TwoPntModelQvec[numpnts(TwoPntModelQvec)-2] )
     344//              InvarData=areaXY(OriginalQvector, OriginalIntensity, Qmin, TwoPntModelQvec[numpnts(TwoPntModelQvec)-2]  )
     345//                      TwoPntModelIntensity*=InvarData/InvarModel
     346//                      CheckDisplayed /W=TwoPhaseSystemData TwoPntModelIntensity
     347//                      if(V_flag==0)
     348//                              AppendToGraph/W=TwoPhaseSystemData  TwoPntModelIntensity vs TwoPntModelQvec
     349//                      endif
     350//                      ModifyGraph lstyle(TwoPntModelIntensity)=9,lsize(TwoPntModelIntensity)=3,rgb(TwoPntModelIntensity)=(2,39321,1)
     351//                      //ModifyGraph mode(TwoPntModelIntensity)=4,marker(TwoPntModelIntensity)=19
     352//                      //ModifyGraph msize(TwoPntModelIntensity)=3
     353//      endif
     354//
     355
     356///******************************************************************************************************************************************
     357///******************************************************************************************************************************************
     358
     359Function IR3T_CalcRadialAveProfile3D(My3DWaveIn)
     360                wave My3DWaveIn
     361                //this calculates radial profile of intensity for 3D Wave
     362                //from the center
     363                //may not be cube... 
     364                variable pcen, qcen, rcen
     365                wavestats/Q My3DWaveIn
     366                pcen = V_maxRowLoc
     367                qcen = V_maxColLoc
     368                rcen = V_maxLayerLoc
     369                variable maxDist=DimSize(My3DWaveIn, 0 )
     370                variable VoxelSize=DimDelta(My3DWaveIn, 0)
     371
     372                //MatrixOp/Free/NTHR=0 My3DTempWv = My3DWaveIn
     373                MatrixOp/Free/NTHR=0 My3DDistacneWv = My3DWaveIn
     374                multithread My3DDistacneWv = sqrt((p-pcen)^2+(r-rcen)^2+(q-qcen)^2)
     375                //make/Free/N=(numpnts(My3DDistacneWv)) My1DValuesWave, My1DDistanceWv
     376                //MatrixOp/Free/NTHR=0 My1DValuesWave = My3DWaveIn                              //this convert 3D wave in 1D wave.
     377                //MatrixOp/Free/NTHR=0 My1DDistanceWv = My3DDistacneWv                  //this convert 3D wave in 1D wave.
     378                //Sort My1DDistanceWv, My1DDistanceWv, My1DValuesWave
     379                make/Free/N=(maxDist) HistogramWvIndx, HistogramWv
     380                Histogram /B={0,1,maxDist}/W=My3DWaveIn/Dest=HistogramWv  My3DDistacneWv
     381                //Wave W_SqrtN          //this does not work with /W flag 
     382                Histogram /B={0,1,maxDist}/Dest=HistogramWvIndx  My3DDistacneWv
     383                HistogramWv /= HistogramWvIndx
     384                Duplicate/O HistogramWv, RadialWaveProfile
     385end
     386               
     387               
    323388///******************************************************************************************************************************************
    324389///******************************************************************************************************************************************
     
    394459//*****************************************************************************************************************
    395460//this function takes as input voxelgram, generated directly or using ConvertToVoxelGram(ListWv, primaryradius)
    396 // VoxelSize here is dimension, which is the size of the voxel... For 3D Solids it is waht model assumed.
     461// VoxelSize here is dimension, which is the size of the voxel... For 3D Solids it is what model assumed.
    397462// If using ConvertToVoxelGram(ListWv, primaryradius) one needs to calculated this from two values:
    398463// Size of the primary particle of the aggregate, close to Rg level 1, may be needs to be slightly modified
     
    400465//   Therefore, if Rg of level 1 is 20 A = assume diameter is 40A, the VoxelSize is actually 4
    401466//   each original "particle" in the Aggregate was replaced with 10 x 10 x 10 voxels, basically we made assumption that diameter of the primary particle is 10.   
    402 Function IR3T_CreatePDF(ThreeDVoxelGram,VoxelSize, NumRSteps, IsoValue, oversample, Qmin, Qmax, NumQSteps)
     467//this function returns PDFIntensityWv and PDFQWv waves with calculated intensity vs Q for the 3D structure.
     468
     469// I M P O R T A N T - verified for spheres on 5-1-2020, seems to work for:
     470// 3D voxelgram, 0 is empty space, 1 is particle
     471// checked on sphere model with voxel sizes 1, 2, and 4 A and it works.
     472// now for 3D solid where we have many particles this does not work. We get again scattering from teh box size and not from particle shapes...
     473// this is probably useful most for Fractal Aggregate and not for 3D solid... 
     474
     475Function IR3T_CreatePDFIntensity(ThreeDVoxelGram, IsoValue, Qmin, Qmax, NumQSteps)
    403476        wave ThreeDVoxelGram
    404         variable VoxelSize, NumRSteps, IsoValue, oversample, Qmin, Qmax, NumQSteps
     477        variable IsoValue, Qmin, Qmax, NumQSteps
    405478        //oversample can be 1 or larger integer (2-4 is sensible). Makes evaluated voxelgram oversample * larger to improve data at high-qs.
    406479        //VoxelSize in [A] per voxel
     
    409482        variable StartTicks=ticks
    410483       
     484        variable VoxelSize = dimdelta(ThreeDVoxelGram,0)                //this is voxelsize from x scaling of 3D wave
     485        variable NumRSteps
     486        variable oversample
     487        //if the 3DVocelgram is sufficiently small, oversample to get better data...
     488        if(dimsize(ThreeDVoxelGram,0)*dimsize(ThreeDVoxelGram,1)*dimsize(ThreeDVoxelGram,2)< 120^3)
     489                oversample = 4
     490        elseif(dimsize(ThreeDVoxelGram,0)*dimsize(ThreeDVoxelGram,1)*dimsize(ThreeDVoxelGram,2)< 220^3)
     491                oversample = 2
     492        else
     493                oversample = 1
     494        endif
    411495        if(oversample>1)                //this helps a lot, takes time, but improves high-q fitting...
    412496                Print "oversampling Data by factor of "+num2str(oversample)
    413497                StartTicks=ticks
    414                 oversample = ceil (oversample)                  //make sure we have integer.
    415                 SetScale/P x 0,1,"", ThreeDVoxelGram            //Interp3D uses wave scaling, so we need to match wave scaling to point numbers...
    416                 SetScale/P y 0,1,"", ThreeDVoxelGram
    417                 SetScale/P z 0,1,"", ThreeDVoxelGram
    418                 make/Free/S/N=(oversample*dimsize(ThreeDVoxelGram,0),oversample*dimsize(ThreeDVoxelGram,1),oversample*dimsize(ThreeDVoxelGram,2)) OversampledThreeDVoxelGram
     498                DUPLICATE/free ThreeDVoxelGram, ThreeDVoxelGramTemp
     499                SetScale/P x 0,1,"", ThreeDVoxelGramTemp                //Interp3D uses wave scaling, so we need to match wave scaling to point numbers...
     500                SetScale/P y 0,1,"", ThreeDVoxelGramTemp
     501                SetScale/P z 0,1,"", ThreeDVoxelGramTemp
     502                make/Free/S/N=(oversample*dimsize(ThreeDVoxelGramTemp,0),oversample*dimsize(ThreeDVoxelGramTemp,1),oversample*dimsize(ThreeDVoxelGramTemp,2)) OversampledThreeDVoxelGram
    419503                //multithread Use3DWave = ThreeDVoxelGram[floor(p/oversample)][floor(q/oversample)][floor(r/oversample)]
    420                 multithread OversampledThreeDVoxelGram = Interp3D(ThreeDVoxelGram, ((p/oversample)), ((q/oversample)),((r/oversample)) )       
     504                multithread OversampledThreeDVoxelGram = Interp3D(ThreeDVoxelGramTemp, ((p/oversample)), ((q/oversample)),((r/oversample)) )   
    421505                MatrixOp/Free/NTHR=0 Use3DWaveThresh = greater(OversampledThreeDVoxelGram,0.55)         //this tresholds the wave so it is smoother than what simple assignment above does...
    422506                wave Use3DWave = Use3DWaveThresh
     
    430514        RadMin = 0.5
    431515        RadMax = ceil(sqrt(3)*max(DimSize(Use3DWave, 0 ), DimSize(Use3DWave, 1),DimSize(Use3DWave, 2 ) ))
    432         Make/O/N=(NumRSteps)/D VoxDistanceWave
    433         VoxDistanceWave = p*(RadMax-RadMin)/NumRSteps   + 1                                             //linear distance space...
     516        NumRSteps = RadMax
     517        Make/free/N=(NumRSteps)/D VoxDistanceWave                                                               //this makes steps in R = 1 A
     518        VoxDistanceWave = p + RadMin                                                                                    //linear distance space...
    434519        //testing shows, that linear radius binning is better. We get better behavior at high q values.
    435520        //IN2G_ConvertTologspacing(VoxDistanceWave, 0.5)                                                        //sets k scaling on log-scale...       
     
    460545        PDFWave[0]=0
    461546        RadiiWave[0]=0
    462        
    463         //this is questionable here - these "radii waves" are actually distances, not radii but steps. Do we need to divide this by two???
    464         RadiiWave/=2
    465         //this surely works very well now for Fractal aggregate, but the hell if I understand  this!!!!!
    466        
     547
    467548        //remove tailing 0 values from calculations. Note: leaves last 0 point in there...
    468549        IR3T_RemoveTailZeroes(PDFwave,RadiiWave)                                                //remove end 0 values, just waste to drag around and plot.
     
    484565        IN2G_ConvertTologspacing(PDFQWv,0)                                                                      //creates log-q spacing in the PDFQWv
    485566        multithread PDFIntensityWv =  IR3T_CalcIntensityPDF(PDFQWv[p],PDFWave,RadiiWave)       
     567        KillWaves/Z PDFWave,RadiiWave
    486568end
    487569//******************************************************************************************************************************************************
     
    538620//******************************************************************************************************************************************************
    539621//******************************************************************************************************************************************************
    540 //threadsafe
    541 threadsafe Function IR3T_CalcIntensityPDF(Qvalue,PDF,Radius)            //AMemiys - my theory powerpoint
    542         variable Qvalue
    543         wave PDF,Radius
    544         Make/Free/N=(numpnts(PDF))/D QRWave
    545         QRWave=sinc(Qvalue*Radius[p])                                                           //(sin(Qvec[p]*Radius))/(Qvec[p]*Radius)               
    546         matrixOP/Nthr=0/Free tempWave = PDF * QRWave
    547         //variable AreaW = areaXY(Radius, TempWave)
    548         return (4*pi*(areaXY(Radius, TempWave)))
    549 end
    550622//******************************************************************************************************************************************************
    551623//******************************************************************************************************************************************************
  • trunk/User Procedures/Irena/IR3_3DTwoPhaseSolid.ipf

    r893 r934  
    1818        //John A Quintanilla, Jordan T Chen, Richard F Reidy and Andrew J Allen Modelling Simul. Mater. Sci. Eng. 15 (2007) S337–S351, doi:10.1088/0965-0393/15/4/S02
    1919
    20 constant useSAXSMorphCode = 0
     20//constant useSAXSMorphCode = 0
    2121
    2222//******************************************************************************************************************************************************
     
    106106        SetVariable TotalNumberOfVoxels,pos={230,290},size={200,16},title="No of Voxels ",noproc, help={"How many voxels is in box, impacts speed!"}
    107107
     108        CheckBox useSAXSMorphCode,pos={100,310},size={250,14},proc=IR3T_InputPanelCheckboxProc,title="USAXS SAXSMorph code (slower)"
     109        CheckBox useSAXSMorphCode,variable=root:Packages:TwoPhaseSolidModel:useSAXSMorphCode, help={"To use GRF generation using SAXSMorph method"}
     110
     111
    108112        //Dist Tabs definition
    109         TabControl TwoPhaseModelTabs,pos={5,320},size={370,280},proc=IR3T_TwoPhaseTabProc
     113        TabControl TwoPhaseModelTabs,pos={5,330},size={370,280},proc=IR3T_TwoPhaseTabProc
    110114        TabControl TwoPhaseModelTabs,tabLabel(0)="1. Extrapolate ",tabLabel(1)="2. Advanced Pars ",tabLabel(2)="3. Results "
    111115
     
    406410
    407411Function IR3T_Calc1DSASData()
    408 
    409                 IR3T_Calc3DTwoPntCorrelation()
    410                 IR3T_CalcTwoPntsCorFIntensity()
     412               
     413                wave My3DWv = root:Packages:TwoPhaseSolidModel:TwoPhaseSolidMatrix
     414                //note, assumes My3DWave has proper x scaling to indicate voxelsize in A.
     415                //this does not work for 3D solid, since it gives scattering from external box size...
     416                //use this for Fractal aggregate and not for 3D solid...
     417                //IR3T_CreatePDFIntensity(My3DWv, 0.5, 0.001, 0.5, 200)
     418                //these are PDF distance calculation results
     419                //Wave PDFIntensityWv
     420                //wave PDFQWv
     421                IR3T_CalcAutoCorelIntensity(My3DWv, 0.5, 0.001, 0.5, 200)
     422                //these are autocorrelation calculated intensities...
     423                Wave AutoCorIntensityWv
     424                Wave AutoCorQWv
     425                Wave OriginalIntensity = root:Packages:TwoPhaseSolidModel:OriginalIntensity
     426                Wave OriginalQvector = root:Packages:TwoPhaseSolidModel:OriginalQvector
     427                variable IntgInt=areaXY(OriginalQvector, OriginalIntensity)
     428                variable ModelIntgInt=areaXY(AutoCorQWv, AutoCorIntensityWv)
     429                AutoCorIntensityWv*=IntgInt/ModelIntgInt
     430               
     431                DOWIndow TwoPhaseSystemData
     432                if(V_Flag)
     433                        DoWIndow/F TwoPhaseSystemData
     434                        //CheckDisplayed /W=TwoPhaseSystemData PDFIntensityWv
     435                        //if(V_flag==0)
     436                        //      AppendToGraph/W=TwoPhaseSystemData/R  PDFIntensityWv vs PDFQWv
     437                        //endif
     438                        CheckDisplayed /W=TwoPhaseSystemData AutoCorIntensityWv
     439                        if(V_flag==0)
     440                                AppendToGraph/W=TwoPhaseSystemData/R  AutoCorIntensityWv vs AutoCorQWv
     441                        endif
     442                        //ModifyGraph lstyle(PDFIntensityWv)=9,lsize(PDFIntensityWv)=3,rgb(PDFIntensityWv)=(1,16019,65535)
     443                        //ModifyGraph mode(PDFIntensityWv)=4,marker(PDFIntensityWv)=19
     444                        //ModifyGraph msize(PDFIntensityWv)=3
     445                        ModifyGraph lsize(AutoCorIntensityWv)=3,rgb(AutoCorIntensityWv)=(3,52428,1)
     446                endif
     447               
    411448end
    412449
     
    809846                        endif
    810847                        //calculate Invariant
    811                         variable Qmax=2         //Qmax=2 hardwired here...
     848                        variable Qmax=5         //Qmax=5 hardwired here...
    812849                        Duplicate/Free ExtrapolatedQvector, ExtrapolatedIntQ2
    813850                        ExtrapolatedIntQ2 = ExtrapolatedIntensity * ExtrapolatedQvector^2
     
    825862                                        Porosity = 0.1
    826863                                        IR3T_SetControlsInPanel()
     864                                else
     865                                        Porosity = (1-sqrt(1-4*Porosity))/2                                     //this is quadratic equation solver
    827866                                endif
    828867                        else
    829868                                ScatteringContrast  = (Invariant / Porosity )*1e-20/(2*pi^2)
    830869                        endif
    831                         Porosity = (1-sqrt(1-4*Porosity))/2                                     //this is quadratic equation solver
    832870                        SurfaceToVolumeRatio = Porosity*(1-porosity)* 1e4*pi*PorodConstant/Invariant    // this is not really S/V, that would be S/V = piB/Q * (phi*(1-phi))
    833871                        if(RKParametersManual==0)               //calculate and set R and K parameters here...
     
    944982        ListOfVariables+="LowQExtrapolationMin;LowQExtrapolationStart;LowQExtrapolationEnd;HighQExtrapolationEnd;HighQExtrapolationStart;HighQExtrapolationMax;"
    945983        ListOfVariables+="PorodConstant;Background;RgValue;VoxelResolution;GizmoFillSolid;"
    946         ListOfVariables+="AchievedVolumeFraction;"
     984        ListOfVariables+="AchievedVolumeFraction;useSAXSMorphCode;"
    947985       
    948986        ListOfStrings="LowQExtrapolationMethod;"
     
    11311169        print "Gamma_alfa(r) calculation time was "+num2str((ticks-startTicks)/60) +" sec"
    11321170        //renormalize
    1133         wavestats/Q PhaseAutocorFnct
     1171        variable MaxValue=wavemax(PhaseAutocorFnct)
    11341172        //Refer to original SAXSMorph java code, and John A. Quintanilla papers, explains why this is normalized as below:
    11351173        // Quintanilla: DebyeAutoCorFnct = S2, and S2(0) = porosity (volume fractio of pores)
     
    11371175        //Java code:      this.gammar[i][1] = ((this.porosity - this.porosity * this.porosity) * gamma_r_val[i] / gamma_r_max + this.porosity * this.porosity);
    11381176        Duplicate/O PhaseAutocorFnct, TwoPointsCorFData                                                                                                                                 //thsi is two points correlation function.
    1139         TwoPointsCorFData = (porosity - porosity^2) * PhaseAutocorFnct[p]/V_max + porosity^2                                    //this converts Two Points Correlation function... .
     1177        TwoPointsCorFData = (porosity - porosity^2) * PhaseAutocorFnct[p]/MaxValue + porosity^2                                 //this converts Two Points Correlation function... .
    11401178        // Makes sense, porosity^2 is random probability that two random points will be both in pore/solid (whatever the minority phase is).
    11411179        //       >>>>>>      this now works, 2-24-2019, but still seem to get different values than SAXSMorph
    11421180
    1143         display PhaseAutocorFnct vs Radii as "PhaseAutocorFnct"
    1144         display TwoPointsCorFData vs Radii as "TwoPointsCorFData"
     1181        //display/K=1 PhaseAutocorFnct vs Radii as "PhaseAutocorFnct"
     1182        //display/K=1 TwoPointsCorFData vs Radii as "TwoPointsCorFData"
    11451183
    11461184                                //**************************************************************
     
    11671205        print "Calculating gR(r) - the two-point correlation function g(r)"             
    11681206        Duplicate/O Radii,AutoCorfnctGr
     1207        NVAR useSAXSMorphCode=root:Packages:TwoPhaseSolidModel:useSAXSMorphCode
    11691208        if(useSAXSMorphCode)
    11701209                print "Calculating g(r) using SAXSMorph code"
     
    11801219                print "Calculating g(r) using code in Quantanilla paper"
    11811220                multithread AutoCorfnctGr = IR3T_ProperCalcgr(alfaValueQ, XiFunctionQuint[p])                           //this is correct way of calculating gr
    1182                 AutoCorfnctGr[0]=1                                                                                                                                                              //first point is 1 by definition and code gets NaN
    1183                 AutoCorfnctGr = numtype(AutoCorfnctGr[p]) == 0 ? AutoCorfnctGr[p] : 1                                                                                   //for log scaled data we have number of nans at the begginign due to really small values, need to be 1 for all such values.
    1184         endif
    1185 
    1186         display AutoCorfnctGr vs Radii as "AutoCorfnctGr"
     1221                AutoCorfnctGr[0]=1                                                                                                                                                                                      //first point is 1 by definition and code gets NaN
     1222                AutoCorfnctGr = numtype(AutoCorfnctGr[p]) == 0 ? AutoCorfnctGr[p] : 1                                           //for log scaled data we have number of nans at the begginign due to really small values, need to be 1 for all such values.
     1223        endif
     1224
     1225        //display/K=1 AutoCorfnctGr vs Radii as "AutoCorfnctGr"
    11871226
    11881227        print "gamma(r) calculation time was "+num2str((ticks-startTicks)/60) +" sec"
     
    12011240        // spectral function is ridiculously noisy, lets smooth it.
    12021241
    1203         display SpectralFk vs Kvalues as "SpectralFk"
     1242        //display/K=1 SpectralFk vs Kvalues as "SpectralFk"
    12041243
    12051244        print "Spectral function calculation time was "+num2str((ticks-startTicks)/60) +" sec"
     
    12151254        //this gets some differeces between SAXSMorph and Igor code, but it may be just different rounding and methods. Not sure which one is actually right here anyway.
    12161255
    1217         display PositiveSpectralFk vs Kvalues as "PositiveSpectralFk"
    1218 
    1219                                 //********************************************
    1220                                 // >>>>   Step 6 from Quintanilla   <<<<<<
    1221                                 //********************************************
    1222                                 //Evaluate G(r), the optimal positive-definite approximation to g(r).
    1223         duplicate/O Radii, PositiveAutoCorfnctGr
    1224         multithread PositiveAutoCorfnctGr =     IR3T_ConvertSpectrFtoGr(Radii[p],PositiveSpectralFk,Kvalues)                   
     1256        //display/K=1 PositiveSpectralFk vs Kvalues as "PositiveSpectralFk"
    12251257                               
    1226         display PositiveAutoCorfnctGr vs Radii as "PositiveAutoCorfnctGr"
     1258        ///display/K=1 PositiveAutoCorfnctGr vs Radii as "PositiveAutoCorfnctGr"
    12271259        startTicks=ticks
    1228         if(!useSAXSMorphCode)
    1229                 print "Using FFT to genrate 2 phase solid, this is not optimized yet"
    1230                 //need to first create
    1231                 //here we can take over by use of FFT...
    1232                 IR3T_UseFFTtoGenerateMatrix(PositiveAutoCorfnctGr,Radii, alfaValue, BoxResolution, BoxSideSize)
    1233         else
     1260        if(useSAXSMorphCode)
    12341261                print "Using SAXSMorph code... "
    12351262                //this is SAXSMorph
     
    12401267                //this is SAXSMoprh way
    12411268                IR3T_GenerateMatrix(Kvalues, SpectralFk, BoxSideSize, BoxResolution, alfaValue)
    1242                 //and this is using FFT
    1243                 //you need FFTGRF.ipf for following step. Its is much faster but the sizes generated are smaller and make not much sense for now...
    1244                 //may be need better - properly fft's gR (bot spectralfk, which has same absolute values and shape, but seem shfited by Rmin...
    1245                 //and ????
    1246                 //      IR3T_GenerateMatrixFFT(SpectralFk, Kvalues, BoxSideSize, BoxResolution, porosity)
     1269        else
     1270                //and this is using FFT based on Quintanilla
     1271                //********************************************
     1272                // >>>>   Step 6 from Quintanilla   <<<<<<
     1273                //********************************************
     1274                //Evaluate G(r), the optimal positive-definite approximation to g(r).
     1275                duplicate/O Radii, PositiveAutoCorfnctGr
     1276                multithread PositiveAutoCorfnctGr =     IR3T_ConvertSpectrFtoGr(Radii[p],PositiveSpectralFk,Kvalues)   
     1277                variable NormalizeVal= wavemax(PositiveAutoCorfnctGr)
     1278                multithread PositiveAutoCorfnctGr = PositiveAutoCorfnctGr/NormalizeVal 
     1279                print "Using FFT to generate 2 phase solid, this is not optimized yet"
     1280                //here we can take over by use of FFT...
     1281                IR3T_UseFFTtoGenerateMatrix(PositiveAutoCorfnctGr,Radii, alfaValue, BoxResolution, BoxSideSize)
    12471282        endif
    12481283        print "GenerateMatrix time is "+num2str((ticks-startTicks)/60) +" sec"
     
    12781313        Wave TwoPhaseSolidMatrix
    12791314        //oversample = BoxResolution>100 ? 2 : 4
     1315        KillWaves/Z RandomField_FFT_Covar_IFFT_norm,AutoCorMatrix,RandomField,CovarFunction3D
     1316
    12801317        print "Done..."
    12811318        resumeUpdate
     
    12961333       
    12971334        //need to add radius 0 = Gr = 1
    1298         Make/O/N=(numpnts(Gr)+1) GRtemp, Radiitemp
    1299         GRtemp[0]=1
    1300         GRtemp[1, ] = Gr[p-1]
     1335        Duplicate/O Gr, GRtemp, Radiitemp
     1336        //GRtemp[0]=1
     1337        GRtemp = Gr/wavemax(Gr)                 //Positive Gr has not been normalized... Need to normalize to 1 at max. 
    13011338        Radiitemp[0]=0
    13021339        Radiitemp[1, ]=Radii[p-1]       
    13031340        //now we need to create Corelation.
    13041341        make/O/N=5000 CovarFunction1D
    1305         variable maxRadgR = radii[numpnts(radii)-1] - 10
    1306         SetScale x -2*maxRadgR,2*maxRadgR,"A", CovarFunction1D
    1307 
     1342        variable maxRadgR = radii[numpnts(radii)-1]
     1343        SetScale x -0.95*maxRadgR,0.95*maxRadgR,"A", CovarFunction1D
     1344
     1345        CovarFunction1D = GRtemp[binarysearchInterp(Radiitemp, abs(x) )]
     1346       
    13081347        //for testing
    1309         Duplicate/O CovarFunction1D, CovarFunction1DGauss, CovarFunction1DOrig
     1348        //Duplicate/O CovarFunction1D, CovarFunction1DGauss, CovarFunction1DOrig
    13101349               
    13111350        //Let's see what using smooth function would do...
    13121351        //Curve fitting with Gaussian works fine, but is not generic enough...
    1313         make/N=4/O W_coef
    1314         W_coef={0,1,0,50}
    1315         K0 = 0;K1 = 1;K2 = 0;K3 = 50;
    1316         CurveFit/G/H="1110" gauss GRtemp /X=Radiitemp   
    1317         CovarFunction1D =  W_coef[0]+W_coef[1]*exp(-((x-W_coef[2])/(W_coef[3]))^2)                              //<<<< is this to convert radius to distance??? I think Covar FUnction needs distance.
     1352        //make/N=4/O W_coef
     1353        //W_coef={0,1,0,50}
     1354        //K0 = 0;K1 = 1;K2 = 0;K3 = 50;
     1355        //CurveFit/G/H="1110" gauss GRtemp /X=Radiitemp
     1356        //CovarFunction1D =  W_coef[0]+W_coef[1]*exp(-((x-W_coef[2])/(W_coef[3]))^2)                            //<<<< is this to convert radius to distance??? I think Covar FUnction needs distance.
    13181357        //end of use of gauss to smooth Gr...
    13191358
    1320                                         //use FFT to remove high osciallations...
     1359                                //use FFT to remove high osciallations...
    13211360                                //      CovarFunction1D = abs(x)<maxRadgR ? GRtemp[binarysearchInterp(Radiitemp, abs(x) )] : 0
    13221361                                //      FFT/OUT=3/DEST=CovarFunction1D_FFT      CovarFunction1D
    13231362                                //      Smooth/EVEN/B 30, CovarFunction1D_FFT
    1324                                         //IFFT/DEST=CovarFunction1D_FFT_SM_IFFT  CovarFunction1D_FFT
     1363                                // IFFT/DEST=CovarFunction1D_FFT_SM_IFFT  CovarFunction1D_FFT
    13251364                                //end, this seems to be a problem...
    13261365               
     
    13361375        //thsi does not fix anything... This is generating too small particles, something is off here.
    13371376        //Need to increase size of the particles twice, but cannot find way to scaqle them up.
    1338 //      SetScale y -1*HalfBox,HalfBox,"A", CovarFunction3D
    1339 //      SetScale x -1*HalfBox,HalfBox,"A", CovarFunction3D
    1340 //      SetScale z -1*HalfBox,HalfBox,"A", CovarFunction3D
    1341 //      SetScale x -1*HalfBox,HalfBox,"A", RandomField
    1342 //      SetScale y -1*HalfBox,HalfBox,"A", RandomField
    1343 //      SetScale z -1*HalfBox,HalfBox,"A", RandomField
    1344 
    1345         multithread CovarFunction3D = CovarFunction1D(sqrt(x^2+y^2+z^2))
     1377        //      SetScale y -1*HalfBox,HalfBox,"A", CovarFunction3D
     1378        //      SetScale x -1*HalfBox,HalfBox,"A", CovarFunction3D
     1379        //      SetScale z -1*HalfBox,HalfBox,"A", CovarFunction3D
     1380        //      SetScale x -1*HalfBox,HalfBox,"A", RandomField
     1381        //      SetScale y -1*HalfBox,HalfBox,"A", RandomField
     1382        //      SetScale z -1*HalfBox,HalfBox,"A", RandomField
     1383
     1384        multithread CovarFunction3D = CovarFunction1D(sqrt(x^2+y^2+z^2))                //this causes error at this time.
    13461385
    13471386        //Now FFT the two components:
    1348         FFT/OUT=1/DEST=RandomField_FFT  RandomField                             // this is random field after FFT, it should be more or less another random Gauss field.
    1349         FFT/OUT=1/DEST=CovarFunction3D_FFT  CovarFunction3D             // this is FFT of covariance function
    1350        
     1387        FFT/OUT=1/Free/DEST=RandomField_FFT  RandomField                                        // this is random field after FFT, it should be more or less another random Gauss field.
     1388        FFT/OUT=1/Free/DEST=CovarFunction3D_FFT  CovarFunction3D                // this is FFT of covariance function
     1389
    13511390        //cannot make this below work somehow....
    1352         //multithread CovarFunction3D_FFT = CovarFunction1D_FFT(sqrt(x^2+y^2+z^2))
    1353 
    1354 
    13551391        matrixOp/NTHR=0/O CovarFunction3D_FFTpwr = powC(CovarFunction3D_FFT, 1/2)       //      << even for 3D this needs to be sqrt, so why is it for 1D not sqrt???
    1356        
    13571392        matrixOp/NTHR=0/O RandomField_FFT_Covar = RandomField_FFT * CovarFunction3D_FFTpwr
    13581393
    13591394        IFFT/DEST=RandomField_FFT_Covar_IFFT  RandomField_FFT_Covar
    1360 
    1361 
     1395        //no normalization lost or needed, when Gr is normalized to 1 at max.
     1396       
    13621397        MatrixOp/NTHR=0/O TwoPhaseSolidMatrix = greater(RandomField_FFT_Covar_IFFT,alfaValue)
     1398
     1399        KillWaves/Z CovarFunction3D_FFTpwr, RandomField_FFT_Covar, RandomField_FFT_Covar_IFFT
    13631400
    13641401        SetScale x 0,BoxSide,"A", TwoPhaseSolidMatrix
     
    14421479Function IR3T_Autocorelate3D(Mat)
    14431480        wave Mat
    1444         FFT/Dest=TmpFFTMat/Free Mat                     //this is fft2(H)
    1445         Duplicate/O/C/Free TmpFFTMat, TmpFFTMatConj, TmpInternal
    1446         multithread TmpFFTMatConj = conj(TmpFFTMat)             //this is conj(fft2(H))
    1447         multithread TmpInternal=TmpFFTMat*TmpFFTMatConj
    1448         IFFT/Dest=AutoCorMatrix TmpInternal
    1449         //AutoCorMatrix/=(dimsize(Mat,0)*dimSize(Mat,1))
    1450         ImageTransform swap3D AutoCorMatrix                     //this is for 3D waves.
    1451         //ImageTransform swap AutoCorMatrix
    1452         wavestats/Q AutoCorMatrix
    1453         AutoCorMatrix = AutoCorMatrix/V_max                             //this scales to max =1 for autocorrelation
     1481        //this autocorelates Mat with itself.
     1482        FFT/Dest=TmpFFTMat/Free Mat                                                                                     //this is slow process... 15 sec in 512^3 wave test
     1483        Duplicate/C/Free TmpFFTMat, TmpFFTMatConj, TmpInternal          //~1 sec in test, MatrixOp is only slightlly faster.
     1484        multithread TmpFFTMatConj = conj(TmpFFTMat)                                             //this is quite fast   
     1485        multithread TmpInternal=TmpFFTMat*TmpFFTMatConj                         //2 sec on test
     1486        IFFT/Dest=AutoCorMatTmp/Free TmpInternal                                                        //this is slow process... 15 sec in 512^3 wave test
     1487        ImageTransform swap3D AutoCorMatTmp                                                             //this is for 3D waves.
     1488        variable MaxValue=WaveMax(AutoCorMatTmp)
     1489        MatrixOp/O/NTHR=0 AutoCorMatrix = AutoCorMatTmp/MaxValue        //this is relatively fast.
     1490        //result is AutoCorMatrix scaled to max=1 in current data folder.
    14541491end
    14551492
     
    14791516        return 4*pi*areaXY(Radius, TempWave)
    14801517end
    1481 
    1482 
     1518///*************************************************************************************************************************************
     1519///*************************************************************************************************************************************
     1520
     1521//threadsafe
     1522threadsafe Function IR3T_CalcIntensityPDF(Qvalue,PDF,Radius)            //AMemiys - my theory powerpoint
     1523        variable Qvalue
     1524        wave PDF,Radius
     1525        Make/Free/N=(numpnts(PDF))/D QRWave
     1526        QRWave=sinc(Qvalue*Radius[p])                                                           //(sin(Qvec[p]*Radius))/(Qvec[p]*Radius)               
     1527        matrixOP/Nthr=0/Free tempWave = PDF * QRWave
     1528        //variable AreaW = areaXY(Radius, TempWave)
     1529        return (4*pi*(areaXY(Radius, TempWave)))
     1530end
    14831531///*************************************************************************************************************************************
    14841532///*************************************************************************************************************************************
     
    16661714        //Kvalues is fk[i][0]
    16671715        //SPectralFK is fk[i][1] in original java code
    1668 
     1716        variable BoxXstep = BoxSideSize/BoxResolution
     1717        variable BoxYstep = BoxSideSize/BoxResolution
     1718        variable BoxZStep = BoxSideSize/BoxResolution
     1719       
    16691720        //description of java code
    16701721        //get Kvalues
     
    16881739   variable IntgNumber = 10000
    16891740        variable xmax, ymax, zmax
    1690         xmax = BoxResolution-1          //num points per side of the box
    1691         ymax = BoxResolution-1
    1692         zmax = BoxResolution-1
     1741        xmax = BoxResolution            //num points per side of the box
     1742        ymax = BoxResolution
     1743        zmax = BoxResolution
    16931744        Make/Free/N=(xmax,3) rmat               //this assumes xmax is at least largest, if not all same???
    16941745        //this is radius matrix, these are dimensions to match with k-vectors randomly generated...
    16951746        // BoxSideSize is real world dimension in Angstroms
    16961747        //
    1697         rmat[][0] = p*BoxSideSize/xmax
    1698         rmat[][1] = p*BoxSideSize/ymax
    1699         rmat[][2] = p*BoxSideSize/zmax
    17001748                //    int xmax = this.BoxResolution;
    17011749                //    int ymax = this.BoxResolution;
     
    17071755                //      rmat[i][2] = (i * 10.0D * this.BoxSideSize / zmax);
    17081756                //    }
     1757        rmat[][0] = 10*p*BoxSideSize/xmax
     1758        rmat[][1] = 10*p*BoxSideSize/ymax
     1759        rmat[][2] = 10*p*BoxSideSize/zmax
    17091760        //in the Java code there is 10.0D* this.BoxSideSize, which is confusing what is it doing...
    17101761        //looks like conversion from A to nm, but the manual indicates it is using A or nm based on Q units and that seems logical. No conversion is done as far as I can say.         
    1711 
     1762        // for some reason, no clue why, this 10 is needed or we get too small features. Confusing, but with 10* it seems to get about the same results as SAXSMorph.
     1763       
    17121764   make/Free/N=3 kvec
    17131765   variable kvecnorm = 0.0
     
    17501802        print "Calculating Gauss random fields, this is the slowest part of the code!"         
    17511803        multithread TwoPhaseSolidMatrix = IR3T_GenGRFUsingCosSaxsMorph(Kn,rmat[p][0],rmat[q][1], rmat[r][2], phin, alpha)       
    1752         SetScale/I x 0,BoxSideSize,"A", TwoPhaseSolidMatrix
    1753         SetScale/I y 0,BoxSideSize,"A", TwoPhaseSolidMatrix
    1754         SetScale/I z 0,BoxSideSize,"A", TwoPhaseSolidMatrix
     1804        SetScale/P x 0,BoxXstep,"A", TwoPhaseSolidMatrix
     1805        SetScale/P y 0,BoxYstep,"A", TwoPhaseSolidMatrix
     1806        SetScale/P z 0,BoxZStep,"A", TwoPhaseSolidMatrix
    17551807
    17561808end
     
    19281980        NVAR GizmoFillSolid =   root:packages:TwoPhaseSolidModel:GizmoFillSolid
    19291981        KillWIndow/Z TwoPhaseSolid3D
    1930 
     1982        variable Xsize, Ysize, ZSize
     1983        Xsize = DimDelta(TwoPhaseSolidMatrix,0)*DimSize(TwoPhaseSolidMatrix,0)
     1984        Ysize = DimDelta(TwoPhaseSolidMatrix,1)*DimSize(TwoPhaseSolidMatrix,1)
     1985        ZSize = DimDelta(TwoPhaseSolidMatrix,2)*DimSize(TwoPhaseSolidMatrix,2)
    19311986        //need to make surfaces here... AGs instructions
    19321987                //2.  I noted that the dimensions are 64x64
     
    20252080        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={10,axisLabel,1}
    20262081        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={11,axisLabel,1}
    2027         ModifyGizmo ModifyObject=axes0,objectType=Axes,property={0,axisLabelText,"300 [A]"}
    2028         ModifyGizmo ModifyObject=axes0,objectType=Axes,property={1,axisLabelText,"300 [A]"}
    2029         ModifyGizmo ModifyObject=axes0,objectType=Axes,property={2,axisLabelText,"300 [A]"}
    2030         ModifyGizmo ModifyObject=axes0,objectType=Axes,property={3,axisLabelText,"300 [A]"}
    2031         ModifyGizmo ModifyObject=axes0,objectType=Axes,property={4,axisLabelText,"300 [A]"}
    2032         ModifyGizmo ModifyObject=axes0,objectType=Axes,property={5,axisLabelText,"300 [A]"}
    2033         ModifyGizmo ModifyObject=axes0,objectType=Axes,property={6,axisLabelText,"300 [A]"}
    2034         ModifyGizmo ModifyObject=axes0,objectType=Axes,property={7,axisLabelText,"300 [A]"}
    2035         ModifyGizmo ModifyObject=axes0,objectType=Axes,property={8,axisLabelText,"300 [A]"}
    2036         ModifyGizmo ModifyObject=axes0,objectType=Axes,property={9,axisLabelText,"300 [A]"}
    2037         ModifyGizmo ModifyObject=axes0,objectType=Axes,property={10,axisLabelText,"300 [A]"}
    2038         ModifyGizmo ModifyObject=axes0,objectType=Axes,property={11,axisLabelText,"300 [A]"}
     2082        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={0,axisLabelText,num2str(Xsize)+" [A]"}
     2083        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={1,axisLabelText,num2str(Xsize)+" [A]"}
     2084        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={2,axisLabelText,num2str(Xsize)+" [A]"}
     2085        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={3,axisLabelText,num2str(Xsize)+" [A]"}
     2086        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={4,axisLabelText,num2str(Ysize)+" [A]"}
     2087        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={5,axisLabelText,num2str(Ysize)+" [A]"}
     2088        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={6,axisLabelText,num2str(Ysize)+" [A]"}
     2089        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={7,axisLabelText,num2str(Ysize)+" [A]"}
     2090        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={8,axisLabelText,num2str(Zsize)+" [A]"}
     2091        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={9,axisLabelText,num2str(Zsize)+" [A]"}
     2092        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={10,axisLabelText,num2str(Zsize)+" [A]"}
     2093        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={11,axisLabelText,num2str(Zsize)+" [A]"}
    20392094        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={0,axisLabelDistance,0}
    20402095        ModifyGizmo ModifyObject=axes0,objectType=Axes,property={1,axisLabelDistance,0}
  • trunk/User Procedures/Irena/IR3_MultiDataPlot.ipf

    r930 r934  
    108108        Button GetHelp,pos={480,10},size={80,15},fColor=(65535,32768,32768), proc=IR3L_ButtonProc,title="Get Help", help={"Open www manual page for this tool"}
    109109        IR3C_MultiAppendControls("Irena:MultiSamplePlot","IR3L_MultiSamplePlotPanel", "IR3L_DoubleClickAction","",0,1)
     110
     111        Button SelectAll,pos={190,680},size={80,15}, proc=IR3L_ButtonProc,title="SelectAll", help={"Select All data in Listbox"}
     112
    110113        //graph controls
    111114        SVAR GraphWindowName=root:Packages:Irena:MultiSamplePlot:GraphWindowName
     
    185188
    186189        SVAR SelectedStyle = root:Packages:Irena:MultiSamplePlot:SelectedStyle
    187         PopupMenu ApplyStyle,pos={260,680},size={400,20},proc=IR3L_PopMenuProc, title="Apply style:",help={"Set tool setting to defined conditions and apply to graph"}
     190        PopupMenu ApplyStyle,pos={280,680},size={400,20},proc=IR3L_PopMenuProc, title="Apply style:",help={"Set tool setting to defined conditions and apply to graph"}
    188191        PopupMenu ApplyStyle,value=#"root:Packages:Irena:MultiSamplePlot:ListOfDefinedStyles",popvalue=SelectedStyle
    189192        Button ApplyPresetFormating,pos={260,710},size={160,20}, proc=IR3L_ButtonProc,title="Apply All Formating", help={"Apply Preset Formating to update graph based on these choices"}
     
    942945                        endif
    943946
     947                        if(stringmatch(ba.ctrlName,"SelectAll"))
     948                                Wave/Z SelectionOfAvailableData = root:Packages:Irena:MultiSamplePlot:SelectionOfAvailableData
     949                                if(WaveExists(SelectionOfAvailableData))
     950                                        SelectionOfAvailableData=1
     951                                endif
     952                        endif
    944953
    945954                        break
Note: See TracChangeset for help on using the changeset viewer.