00001
00002 c
00003 c Copyright (c) 1986,1987,1988,1989,1990,1991,1992,1993,
00004 c by Steve McMillan, Drexel University, Philadelphia, PA.
00005 c
00006 c All rights reserved.
00007 c
00008 c Redistribution and use in source and binary forms are permitted
00009 c provided that the above copyright notice and this paragraph are
00010 c duplicated in all such forms and that any documentation,
00011 c advertising materials, and other materials related to such
00012 c distribution and use acknowledge that the software was developed
00013 c by the author named above.
00014 c
00015 c THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
00016 c IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
00017 c WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
00018 c
00019
00020 program energies
00021 c
00022 parameter (XLEN = 5., YLEN = 5., YMAX = 12.)
00023 parameter (rc = 1., rh = 50., rt = 500.)
00024 c
00025 parameter (N = 17, ALPHA = .4)
00026 dimension e(0:N+1),decm(N)
00027 c
00028 dimension phi(500),rapo(500)
00029 c
00030 e(0) = 0.
00031 e(1) = 1.
00032 do 10 i=2,N+1
00033 10 e(i) = e(i-1)*(1.+ALPHA)
00034 c
00035 do 20 i=1,N
00036 20 decm(i) = (e(i+1) - e(i))/6.
00037 c
00038 ddphi = .05
00039 dphi = -ddphi
00040 nph = 0
00041 c
00042 30 dphi = dphi + ddphi
00043 nph = nph + 1
00044 phi(nph) = dphi
00045 rapo(nph) = log10(apocenter(rc,rh,rt,dphi))
00046 if (rapo(nph).lt.log10(rt).and.nph.lt.500) go to 30
00047 c
00048 c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
00049 c
00050 call mcinit
00051 call nobounds
00052 call sethts(.25,.2)
00053 call plot(4.25,1.,-3)
00054 c
00055 call strpos(1.,1.)
00056 call simbol(XLEN,YLEN+1.2,.15,'Figure 1%%',0.,999)
00057 call clrstr
00058 c
00059 call weight(10)
00060 call eframe(rc,rt,XLEN,-1,'!1.1!r\\a',0.,YMAX,YLEN,2,' ')
00061 c
00062 call simsize(.25,'@D@f',4,dr1,ds)
00063 call simsize(.25,'@s^2',4,dr2,ds)
00064 dr = max(dr1,dr2)
00065 s = 7.5*YLEN/YMAX
00066 call plot(-.25-dr,s,3)
00067 call plot(-.25,s,2)
00068 c
00069 call strpos(1.,0.)
00070 call simbol(-.25,s+.05,.25,'@D@f',0.,4)
00071 call strpos(1.,1.)
00072 call simbol(-.25,s-.05,.25,'@s^2',0.,4)
00073 call clrstr
00074 call mline(rapo,phi,nph,0,0,0.)
00075 c
00076 c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
00077 c
00078 do 60 i=7,N
00079 if (i.eq.8) go to 60
00080 s = decm(i)*YLEN/YMAX
00081 if (s.le.YLEN) then
00082 call dplot(0.,s,3)
00083 do 55 j=1,nph
00084 if (phi(j).gt.decm(i)) then
00085 r = rapo(j-1) + (rapo(j)-rapo(j-1))
00086 & *(decm(i)-phi(j-1))
00087 & /(phi(j)-phi(j-1))
00088 call frinches(r,ss,r,ss)
00089 go to 56
00090 end if
00091 55 continue
00092 56 call setpat(1,2,1,2)
00093 call dplot(r,s,2)
00094 call setpat(1,2,1,2)
00095 call dplot(r,0.,3)
00096 call dplot(r,s,2)
00097 end if
00098 60 continue
00099 c
00100 c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
00101 c
00102 call plot(-2.25,0.,-3)
00103 call strpos(0.,.5)
00104 do 40 i=1,N
00105 s = decm(i)*YLEN/YMAX
00106 if (s.le.YLEN) then
00107 call plot(0.,s,3)
00108 call plot(.6,s,2)
00109 if (i.eq.7.or.i.ge.9) then
00110 call nomber(.65,s,.1,float(i),0.,-1)
00111 call setpat(1,2,1,2)
00112 call dplot(.925,s,3)
00113 call dplot(1.325,s,2)
00114 end if
00115 end if
00116 40 continue
00117 c
00118 call strpos(.5,0.)
00119 call simbol(.6,YLEN,.225,'%E\\C\\M/@s^2%%',0.,999)
00120 c
00121 r1 = 1.5
00122 s1 = decm(7)*YLEN/YMAX
00123 c
00124 c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
00125 c
00126 call plot(-r1,0.,-3)
00127 call strpos(0.,.5)
00128 do 50 i=0,N
00129 s = e(i)*YLEN/YMAX
00130 if (s.le.YLEN) then
00131 call plot(0.,s,3)
00132 call plot(.6,s,2)
00133 if (i.gt.0) then
00134 call nomber(.65,s,.1,float(i),0.,-1)
00135 else
00136 call weight(1)
00137 dr = .6/15
00138 do 45 r = dr,.6,dr
00139 call plot(r-dr,s-2.5*dr,3)
00140 call plot(r,s,2)
00141 45 continue
00142 call weight(10)
00143 end if
00144 end if
00145 50 continue
00146 c
00147 call strpos(.5,0.)
00148 call simbol(.34,YLEN,.2,'E\\b\\i\\n/kT%%',0.,999)
00149 c
00150 r0 = .8
00151 c
00152 ss8 = e(8)*YLEN/YMAX
00153 call plot(r0,ss8,3)
00154 call plot(r0+.05,ss8,2)
00155 ss7 = e(7)*YLEN/YMAX
00156 call plot(r0+.05,ss7,2)
00157 call plot(r0,ss7,2)
00158 c
00159 r0 = r0 + .05
00160 s0 = .5*(ss7+ss8)
00161 r1 = r1 - .05
00162 c
00163 frac = .985
00164 r1 = r0 + frac*(r1-r0)
00165 s1 = s0 + frac*(s1-s0)
00166 c
00167 call plot(r0,s0,3)
00168 call plot(r1,s1,2)
00169 frac = .1
00170 r0 = r1 - frac*(r1 - r0)
00171 s0 = s1 - frac*(s1 - s0)
00172 call arrow(r0,s0,r1,s1,1)
00173 c
00174 c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
00175 c
00176 call mcquit
00177 c
00178 end
00179
