Ising_2D/plot_range.py at master · mtgraves/Ising_2D · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
# =============================================================================
# Script for plotting the data generated from the tempRangeAnalyze.py script.
# This takes reduced file and plots magnetization, energy, and specific heat
# vs. temperature.
#
# Author:           Max Graves
# Last Revision:    21-MAY-2013
# =============================================================================

import subprocess, os, sys, argparse, glob
import DerApproximator as derr
import pylab as pl
import numpy as np
from scipy import interpolate

# =============================================================================
def parseCMD():
    desc= ('Ising 2D script for plotting phase transition.  Takes as its \
            argument the directory holding all of the data files that you ,\
            want to plot.')

    parser = argparse.ArgumentParser(description=desc)
    parser.add_argument('direcN', help='name of reduced file of interest')
    return parser.parse_args()

# =============================================================================
def numericalDer(x,y):
    '''
    calculate dy by 4-point center differencing using array slices

    frac{y[i-2] - 8y[i-1] + 8[i+1] - y[i+2]}{12h}

    y[0] and y[1] must be defined by lower order methods
    and y[-1] and y[-2] must be defined by lower order methods
    '''

    dy = np.zeros(y.shape,np.float) # we know it will be this size
    h = x[1]-x[0] # assumes the points are evenly spaced
    dy[2:-2] = (y[0:-4] - 8*y[1:-3] + 8*y[3:-1] - y[4:])/(12.*h)

    dy[0] = (y[1]-y[0])/(x[1]-x[0])
    dy[1] = (y[2]-y[1])/(x[2]-x[1])
    dy[-2] = (y[-2] - y[-3])/(x[-2] - x[-3])
    dy[-1] = (y[-1] - y[-2])/(x[-1] - x[-2])

    return dy

# =============================================================================
def main():

    carolinaBlue = '#56A0D3'

    # load in reduced data from file
    args = parseCMD()
    fileName = args.direcN
    temps, Es, Ms, Cv = pl.loadtxt(fileName, unpack=True)

    # === take numerical/ analytical derivatives =====
    k = [(temps[i],Es[i],Cv[i]) for i in np.argsort(temps)]

    TsSorted, EsSorted, CvSorted = pl.array([]), pl.array([]), pl.array([])
    for tup in k:
        TsSorted = pl.append(TsSorted, tup[0])
        EsSorted = pl.append(EsSorted, tup[1])
        CvSorted = pl.append(CvSorted, tup[2])

    tempsMore = pl.arange(TsSorted[0],TsSorted[-1], 0.01)
    tck = interpolate.splrep(TsSorted, EsSorted, s=30)
    Esmore = interpolate.splev(tempsMore,tck,der=0)

    dy = numericalDer(TsSorted,Esmore)
    # ================================================

    # plot energy vs. temp
    fig1 = pl.figure(1)
    p1 = fig1.add_subplot(111)
    pl.plot(temps,Es, marker='o', linewidth=0, ms=6,
            markerfacecolor=carolinaBlue, markeredgecolor='Black')
    pl.grid(True)
    pl.xlabel('Temperature '+r'$[K]$', size=20)
    pl.ylabel('Energy', size=20)

    # plot Cv vs. temp using fluctuations
    fig2 = pl.figure(2)
    p2 = fig2.add_subplot(111)
    pl.plot(temps,Cv, marker='o', linewidth=0,
            markerfacecolor='Lime', markeredgecolor='Black')
    pl.grid(True)
    pl.xlabel('Temperature '+r'$[K]$', size=20)
    pl.ylabel('Specific Heat (using fluctuations)', size=20)

    # plot magnetization vs. temp
    fig3 = pl.figure(3)
    p3 = fig3.add_subplot(111)
    pl.plot(temps,pl.absolute(Ms), marker='o', linewidth=0,
            markerfacecolor='Cyan', markeredgecolor='Black')
    pl.grid(True)
    pl.xlabel('Temperature '+r'$[K]$', size=20)
    pl.ylabel('abs(Magnetization)', size=20)

    # plot Cv vs. temp using fluctuations
    fig4 = pl.figure(4)
    p4 = fig4.add_subplot(111)
    pl.plot(tempsMore,dy)
    pl.grid(True)
    pl.xlabel('Temperature '+r'$[K]$', size=20)
    pl.ylabel('Specific Heat (spline fitting)', size=20)

    pl.show()

# =============================================================================
if __name__=='__main__':
    main()