(* file xtalmr/matrix.math *) (* mathematica file -- Interprets the 3-D matrices written by *) (* XREFIn SEARch ROTAtion and XREFin SEARch TRANslation *) Clear[CP]; Off[General::spell1]; \$DefaultFont={"Times-Roman",13}; (*===>*) (* Specify the file name of the matrix written by X-PLOR.*) << rotation.3dmatrix ; (* With the following statement, one can permute the indices of the *) (* 3d-matrix. This is useful if one wants to plot sections of the *) (* first, second, or third level. Example: *) (* perm=(2,1,3); means that one is going to plot 2-sections, where *) (* each 2-section has the 1-coordinate along the x-axis and *) (* the 3-coordinate along the y-axis. *) (* The meaning of the first, second, and third level depends on the *) (* particular application; for example, in the case of spherical *) (* polar angles in XREFin SEARch ROTAtion, it is phi, *) (* kappa, psi. In this case, one might want to use perm={2,1,3}; *) (* one must not permute the indices for matrices with variable *) (* second and third dimensions. This is the case for Lattman's *) (* angle sampling in XREFin SEARch ROTAtion. *) (*===>*) perm={1,2,3}; If[ perm != {1,2,3}, {iperm=perm;iperm[[perm[[1]]]]=1; iperm[[perm[[2]]]]=2;iperm[[perm[[3]]]]=3; rf=Transpose[rf,iperm]; min=min[[perm]];max=max[[perm]];var=var[[perm]];}, {}]; nz=Dimensions[rf][[1]]; (* In the following, one loops through all sections as a function *) (* of the first level. This can be pretty verbose. To plot a *) (* particular section near "z", one can compute the index ii by *) (* using the following statement: *) (* ii=Round[(z-min[[1]])/(max[[1]]-min[[1]])*(nz-1)+1]. *) (* Of course, "z" should be set before issuing this statement. *) (* Then remove the Do loop and plot the single section. *) (*===>*) Do[ { section=Transpose[rf[[ii]]]; (* Mathematica needs a square matrix to make contour plots.*) (* One has to pad the rectangular matrices. *) {ny,nx}=Dimensions[section]; nmax=Max[nx,ny]; If[ nz == 1, {label=StringJoin[var[[1]],"=",ToString[min[[1]]]] }, {label=StringJoin[var[[1]],"=",ToString[(ii-1)/(nz-1) (max[[1]]-min[[1]])+min[[1]] ] ]; } ]; Which[ ny==1, (* One-x-dimensional plot.*) { column0= Table[min[[2]]+(i-1)(max[[2]]-min[[2]])/nx, {i,1,nx}]; xyp=Transpose[{column0,(Flatten[section]-rave)/rsigma}]; ListPlot[xyp, PlotLabel->label, PlotRange->{ {min[[2]],max[[2]]},{ 0 ,(rmax-rave)/rsigma} }, AxesLabel->{ var[[2]], "sigma" }, PlotJoined->True, (*===>*) Ticks->{Range[min[[2]],max[[2]],50.], Range[0,(rmax-rave)/rsigma+0.5,1.]}, AxesStyle->{PostScript ["/Times-Roman findfont 13 scalefont setfont"]}]; }, nx==1, (* One-y-dimensional plot.*) { column0= Table[min[[3]]+(i-1)(max[[3]]-min[[3]])/ny, {i,1,ny}]; xyp=Transpose[{column0,(Flatten[section]-rave)/rsigma}]; ListPlot[xyp, PlotLabel->label, PlotRange->{ {min[[3]],max[[3]]},{ 0 ,(rmax-rave)/rsigma} }, AxesLabel->{ var[[3]], "sigma" }, PlotJoined->True, (*===>*) Ticks->{Range[min[[3]],max[[3]],50.], Range[0,(rmax-rave)/rsigma+0.5,1.]}, AxesStyle->{PostScript ["/Times-Roman findfont 13 scalefont setfont"]}]; }, True,{ (* Two-dimensional plot.*) If[ nmax > ny, { xxx=Table[rave, {j, 1, nmax nmax -nx ny}]; section=Partition[Flatten[Append[section, xxx]], nmax]; pymax=nmax/ny ( max[[3]] - min[[3]] ) + min[[3]]; pxmax=max[[2]]; } , { xxx=Table[rave, {j, 1, nmax nmax -nx ny}]; section=Transpose[section]; section=Partition[Flatten[Append[section, xxx]], nmax]; section=Transpose[section]; pymax=max[[3]]; pxmax=nmax/nx ( max[[2]] - min[[2]] ) + min[[2]]; } ]; CP[ii]=ListContourPlot[section, MeshRange->{{min[[2]],pxmax},{min[[3]],pymax}}, PlotLabel->label, (*===>*) (* Specify plotrange: levels starting at *) (* 2 sigma above the mean. *) PlotRange->{ rave + 2 rsigma , rmax }, Axes->True, AxesLabel->{ var[[2]], var[[3]] }, (*===>*) Ticks->{Range[min[[2]],max[[2]],50.], Range[min[[3]],max[[3]],50.]}, (*===>*) (* Specify number of contour levels.*) Contours->10, ContourShading->False, ContourSmoothing->None]; }]; },{ii,nz}];