The package phasePlot3D

© 2000 Heiko Feldmann

The package phasePlot3D serves to create three-dimensional plots of a function with discontinuities (e.g phase boundaries).
This may be for example some thermodynamic parameter plotted in the phase diagram,
but also any function with branch cuts in the complex plane.

Introduction

First, one needs to define one function which gives information about the "phase diagram":
which function to plot in which region. It is typically called phaseSelection[{x,y}] and 
returns integer values 1,2,3...n, where n is the number of different "phases". Typically,
it will consist of an "If"-command or a "Which"-command.
An example is shown here, defining two circles around {0,0} and {2,0}
as two phases and the rest as the third phase:

\(phaseSelection[{x_, y_}] := 
      Which[
		{x, y} . {x, y} < 0.8, 1, 

        		{x - 2, y} . {x - 2, y} < 0.8, 2, 
		True, 3]\)

The second ingredient are the functions to be plotted in the different regions. I called them phase1[{x,y}] etc.
As usual, one has to define the ranges of the variables in x- and y- direction.
Finally, phasePlot3D uses the same options as Graphics3D.
In total, the complete call has the form

phasePlot3D[
      phaseSelection, {phase1, phase2 ,  ... }, {x, xmin, xmax}, {y, ymin, 
        ymax}]

Two things should be noticed:
1, The functions are used as "pure functions", i.e. without explicitely giving their arguments
2, The functions are expected to be of the form function[{x,y}], with the position {x,y} as a list of length 2 as the only argument.
Apart from this, phasePlot3D can be used as an intuitive extention of Mathematica.

Example 1: What dithering can do to a cylinder

Needs["Graphics`phasePlot3D`", "phasePlot3D.m"];

\(phase1[{x_, y_}] := 1; 
phase2[{x_, y_}] := \(-1\); 

    phase3[{x_, y_}] := 0; 

    phaseSelection[{x_, y_}] := 
      Which[
		{x, y} . {x, y} < 0.8, 1, 

        		{x - 2, y} . {x - 2, y} < 0.8, 2, 
		True, 3]; 

    phasePlot3D[
      phaseSelection, {phase1, phase2, phase3}, {x, \(-1\), 3}, {y, \(-1\), 
        1}, BoxRatios\  -> \ Automatic]; \)

This would be the same function as a standard three-dimensional plot:

Plot3D[
      Switch[phaseSelection[{x, y}], 1, phase1[{x, y}], 2, phase2[{x, y}], 3, 
        phase3[{x, y}]], {x, -1, 3}, {y, -1, 1}, 
      BoxRatios  ->  Automatic];

Example 2: The tricritical phase diagram of the BEG-model

The BEG-model is a spin-1 Ising model which simulates the thermodynamic behavior of (He^3) - (He^4) mixtures at low temperatures.
It was one of the first systems discovered that have a tricritical point. Details can be found in
M. Blume, V. J. Emery and R. B. Griffiths, Phys. Rev. A 4, p.1071
The external parameters are T=Temperature and Δ=chemical potential.
Shown in the (T,Δ).plane is the concentration of (He^3), which serves as an order parameter.

The equations and a tiny bit of physics

To determine the first order line, we need the free energy f and the magnetisation m.
The transition line is obtained by equating the two solutions with M=0 and M≠0 below the tricritical point, while for the second order transition line there exists an analytic formula.

Needs["Graphics`phasePlot3D`", "phasePlot3D.m"]

\(m[T_, Δ_] := \ 
      M\  /. FindRoot[
          M == 2\ Sinh[M/T]\/\(Exp[Δ/T] + 2  Cosh[M/T]\), {M, 
            1}]; 
f[T_, Δ_, M_] := 
      T\ Log[Exp[Δ/T] + 2  Cosh[M/T]] - M\^2/2; 

    firstOrderPoints = 
      Table[{T, 
          Δ\  /. 
            FindRoot[
              f[T, Δ, m[T, Δ]] - 
                f[T, Δ, 0], {Δ, 0.4, 0.46}]}, {T, 
          1. /60, 1. /3 - 1. /60, 1. /60}]; 

    firstOrderLine\  = \ 
      Interpolation[
        Join[firstOrderPoints, {{0. , 0.5}, {1. /3, 2/3\ Log[2]}}]]; 

    secondOrderLine[T_] := T\ Log[2 \((1 - T)\)/T]; 

    transition[T_] := 
      If[T\  < \ 1. /3, firstOrderLine[T], secondOrderLine[T]]; \)

This function determines the thermodynamic phase:

phaseSelection[{Δ_, T_}] := 
      If[Δ > transition[T], 1, 2];

The two formulas for the (He^3)-concentration x,
for M=0 (phase1) and for M≠0 (phase2):

\(phase1[{Δ_, T_}] := 
      1 - 2\/\(Exp[Δ/T] + 2\); 

    phase2[{Δ_, T_}] := 
      1 - m[T, Δ]\ Coth[m[T, Δ]/T]; \)

The 3D plot

phasePlot3D[
      phaseSelection, {phase1, phase2}, {x, 0.1, 0.7}, {y, 0.1, 0.7}, 
      ViewPoint -> {-3,  1,  2}];

Again, one might compare it with the standard plot:

Plot3D[
      If[phaseSelection[{x, y}] == 1, phase1[{x, y}], phase2[{x, y}]], {x, 
        0.7, 0.1}, {y, 0.1, 0.7}, ViewPoint -> {-3,  1,  2}];