User:Lalluviamola/phase space graphs

From Wikimedia Commons, the free media repository
Jump to navigation Jump to search
Phase space of cubic potential

This is the code for doing this graph.

width=12
height=18.5
size width height

!Define height of the graphs
graph_height = height/2

include "graphutil.gle"
set font texcmr

! Hamiltonian of the unidimensional system
! H(q,p)=p^2/(2m) + q/a*( 1+ (q/b)^2 ) 
! H(q,p)=E=T(p)+U(q)


! parameters: different values should no change anything
a=1
b=1
! mass of the particle
m=1
! q where u(q) is an extrem  
q_min_rel = -b/sqrt(3)
q_max_rel = b/sqrt(3)
! u extrems
u_min_rel = -2*b/(a*3*sqrt(3))
u_max_rel = -u_min_rel

! Set extrems of the graphic
q_min = -1.3*b!1.5*q_min_rel
q_max = -1.1*q_min
u_min = -0.7*b/a
u_max = -u_min
p_min = 3.5*u_min_rel
p_max = -p_min

! energy references
e1=1.2*u_min_rel
e2=u_min_rel
e3=0.4*u_min_rel
e4=0.7*u_max_rel
e5=u_max_rel
e6=1.2*u_max_rel

! Some particular values
q3=-0.150548*b
u3=q3/a*( 1+ (q3/b)^2)

q4=0.253198*b
u4=q4/a*( 1+ (q4/b)^2)

amove 0 graph_height

begin graph

   size width graph_height
   math
   title "U(q)=q/a+q^3/b^2"
   ylabels off
   yticks off
   ysubticks off 
   xaxis min q_min max q_max dticks b nticks 3
   xnames "-b" "0" "b"
   yaxis min u_min max u_max

   let d1 = x/a-(x^3)/(a*b^2)
   ! Draw the line
   d1 line color red

end graph

!Draw q-references
!set lstyle 4 just lc
!graph_vline q_min_rel
!graph_vline q_max_rel
!graph_text q_min_rel u_min_rel label "-2b/3(3)^{1/2}a" dx -3

!Draw E-references
set lstyle 1 just lc
set color seagreen
graph_hline e1
amove xg(q_min) yg(e1+u_max/20)
tex "$E_1$"

set color green 
graph_hline e2
amove xg(q_min) yg(e2+u_max/20)
tex "$E_2$"

set color steelblue
graph_hline e3
amove xg(q_min) yg(e3+u_max/20)
tex "$E_3$"

set color deeppink
graph_hline e4
amove xg(q_min) yg(e4+u_max/20)
tex "$E_4$"

set color indigo
graph_hline e5
amove xg(q_min+q_max/8) yg(e5+u_max/20)
tex "$E_5$"

set color chocolate
graph_hline e6
amove xg(q_min+q_max/8) yg(e6+u_max/20)
tex "$E_6$"

!Reset color
set color black

! Grab some absolute positions
y_min_graph_top=yg(u_min_rel)
y_max_graph_top=yg(u_max_rel)
x_min_graph_top=xg(q_min_rel) 
x_max_graph_top=xg(q_max_rel)

q3_top = xg(q3)
u3_top = yg(u3)

q4_top = xg(q4)
u4_top = yg(u4)

! Draw name of the x variable 
amove xg(q_max+b/10) yg(0-b/40)
text q

! Go to the second graph
amove 0 0


begin graph

    ! Define the curves
    E=e1
    let d2 =  sqrt(2*m)*sqrt( E-(x/a)*(1-(x/b)^2) )
    let d3 = -sqrt(2*m)*sqrt( E-(x/a)*(1-(x/b)^2) )
    E=e2
    let d4 =  sqrt(2*m)*sqrt( E-(x/a)*(1-(x/b)^2) )
    let d5 = -sqrt(2*m)*sqrt( E-(x/a)*(1-(x/b)^2) )
    E=e3
    let d6 =  sqrt(2*m)*sqrt( E-(x/a)*(1-(x/b)^2) )
    let d7 = -sqrt(2*m)*sqrt( E-(x/a)*(1-(x/b)^2) )
    E=0
    let d8 =  sqrt(2*m)*sqrt( E-(x/a)*(1-(x/b)^2) )
    let d9 = -sqrt(2*m)*sqrt( E-(x/a)*(1-(x/b)^2) )
    E=e4
    let d10 =  sqrt(2*m)*sqrt( E-(x/a)*(1-(x/b)^2) )
    let d11 = -sqrt(2*m)*sqrt( E-(x/a)*(1-(x/b)^2) )
    E=e5
    let d12 =  sqrt(2*m)*sqrt( E-(x/a)*(1-(x/b)^2) )
    let d13 = -sqrt(2*m)*sqrt( E-(x/a)*(1-(x/b)^2) )
    E=e6
    let d14 =  sqrt(2*m)*sqrt( E-(x/a)*(1-(x/b)^2) )
    let d15 = -sqrt(2*m)*sqrt( E-(x/a)*(1-(x/b)^2) )

   ! Define the graph
   size width graph_height
   math
   ylabels off
   yticks off
   ysubticks off 
   title "p(q)"
   xaxis min q_min max q_max dticks b

    d2  line color seagreen
    d3  line color seagreen
    d4  line color green
    d5  line color green
    d6  line color steelblue
    d7  line color steelblue
    d8  line color black
    d9  line color black
    d10 line color deeppink
    d11 line color deeppink
    d12 line color indigo
    d13 line color indigo
    d14 line color chocolate
    d15 line color chocolate

   xnames "-b" "0" "b"
   yaxis min p_min max p_max 


end graph

! Draw some arrows for showing how system evolves
set lstyle 9 color black

!Central arrow
amove xg(q_min_rel) yg(0)
arc b/2 40 340 arrow start

!Above left arrow
amove xg(q_min_rel) yg(1.15*b)
arc b 40 130 arrow start

!Below right arrow
amove xg(q_max_rel) yg(-b)
arc b 60 130 arrow end

!Right arrow
amove xg(2*b) yg(0)
arc 2*b 160 200 arrow start

! Draw name of the x variable 
amove xg(q_max+b/10) yg(0-b/40)
set color black
text q

! Grab absolute positions
y_min_graph_down=yg(0)
y_max_graph_down=yg(0)
x_min_graph_down=x_min_graph_top
x_max_graph_down=x_max_graph_top

q3_down=xg(q3)
q4_down=xg(q4)

! Draw lines joining both graphs
set lstyle 4 just lc

amove x_min_graph_top   y_min_graph_top
aline x_min_graph_down y_min_graph_down

amove x_max_graph_top   y_max_graph_top
aline x_max_graph_down y_max_graph_down

set color steelblue
amove q3_top u3_top
aline q3_down y_max_graph_down

set color deeppink
amove q4_top u4_top
aline q4_down y_max_graph_down

! Draw a cross for refering the stable equilibrium
set lstyle 1
amove xg(q_min_rel) yg(0)
set color green
marker asterisk 

! Draw a cross for refering the unstable equilibrium
amove xg(q_max_rel) yg(0)
set color indigo
marker asterisk