(* 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}];