博文

关于画图软件Supermongo的一点总结

已有 8701 次阅读 2010-5-6 16:20 |个人分类:总结|系统分类:科研笔记| 编程, 画图

_{Supermongo是一些天文工作者常用的画图工具，曾经可以免费获取，但现在好像要收费了。我原来下载过一个版本，现在也不时用一下。

Supermongo的安装可以参考安装文件夹里的sm.install文件。基本就是
[xx@xx]$ set_opts
这一步要选一些参数，设定安装路径（包括可被C和Fortran调用的库的路径）。}[xx@xx]$ make
[xx@xx]$ make install
至此就安装完毕了。具体的使用没有太多可说的，用C和Fortran调用Supermongo的库有时候对计算有很大帮助，这样一来就可以在计算的同时可视化了，也就是可以监控计算的过程，出错了可以及时发现。我经常求解微分方程，这个功能对我比较有用。

下面就放一段程序作为示例。编译的时候要做如下操作（路径根据库的位置更改）
[xx@xx]$ g77 -c bs.f
[xx@xx]$ g77 -o g77 bs.o -L/home/lqian/local/lib/ -lplotsub -ldevices -lutils -L/usr/X11R6/lib/ -lX11 -o bs
然后就生成了可执行文件bs。

PROGRAM bs
INTEGER nvar,nok,nbad
REAL*8 EPS
PARAMETER (EPS=1.d-12,nvar=4)
REAL*8 hmin,h1,x1,x2,ystart(nvar)
REAL smx,smmu,smnu,smsigma,smsigmap
REAL*8 Omega,Lambda
EXTERNAL bsstep,drv
COMMON/para/ Omega, Lambda

Omega=2.0d0
Lambda=0.0d0

cccccccccccc for plot on screen:
if(sm_device("x11 -bg 0").lt.0)then
cccccccccccc for plot in file test.ps:
c if(sm_device("postfile :SY@: :OF@: equip.ps ").
c 1 lt.0)then
print*,'Cannot open required device'
stop
endif

call sm_graphics
call sm_defvar("TeX_strings","1")
call sm_expand(1.2)

call sm_limits(0.0,100.0,-0.1,0.1) ! set x,y axises

call sm_location(3500,32000,3500,32000)
call sm_box(1,2,0,0)
call sm_xlabel('x')
call sm_ylabel('...')

call sm_alpha

mu0=0.0d0
nu0=0.0d0
sigma0=0.1d0
sigmap0=0.0d0

x1=0.0D0
x2=0.01D0

ystart(1)=mu0
ystart(2)=nu0
ystart(3)=sigma0
ystart(4)=sigmap0

hmin=0.0D0

OPEN(1,FILE='bs.txt')
DO 5 WHILE(x2.LT.1000.0D0)
h1=(x2-x1)/10000.0

CALL odeint(ystart,nvar,x1,x2,EPS,h1,hmin,nok,nbad,drv,bsstep)
PRINT *,x2,ystart(1),ystart(2),ystart(3),ystart(4)
c WRITE(1,*)ystart(1),ystart(2)

smx=x2
smmu=ystart(1)
smnu=ystart(2)
smsigma=ystart(3)
smsigmap=ystart(4)
call sm_ctype('blue')
call sm_ptype(41.,0)
call sm_points(smx,smmu,1)
call sm_ctype('red')
call sm_ptype(41.,0)
call sm_points(smx,smnu,1)
call sm_ctype('yellow')
call sm_ptype(41.,0)
call sm_points(smx,smsigma,1)
call sm_ctype('green')
call sm_ptype(41.,0)
call sm_points(smx,smsigmap,1)

x1=x2
x2=x1+0.01
5 CONTINUE
CLOSE(1)

PAUSE
END

SUBROUTINE odeint(ystart,nvar,x1,x2,eps,h1,hmin,nok,nbad,derivs,
*rkqs)
INTEGER nbad,nok,nvar,KMAXX,MAXSTP,NMAX
REAL*8 eps,h1,hmin,x1,x2,ystart(nvar),TINY
EXTERNAL derivs,rkqs
PARAMETER (MAXSTP=100000,NMAX=50,KMAXX=200,TINY=1.e-30)
INTEGER i,kmax,kount,nstp
REAL*8 dxsav,h,hdid,hnext,x,xsav,dydx(NMAX),xp(KMAXX),y(NMAX),
*yp(NMAX,KMAXX),yscal(NMAX)
COMMON /path/ kmax,kount,dxsav,xp,yp
x=x1
h=sign(h1,x2-x1)
nok=0
nbad=0
kount=0
do 11 i=1,nvar
y(i)=ystart(i)
11 continue
if (kmax.gt.0) xsav=x-2.*dxsav
do 16 nstp=1,MAXSTP
call derivs(x,y,dydx)
do 12 i=1,nvar
yscal(i)=abs(y(i))+abs(h*dydx(i))+TINY
12 continue
if(kmax.gt.0)then
if(abs(x-xsav).gt.abs(dxsav)) then
if(kount.lt.kmax-1)then
kount=kount+1
xp(kount)=x
do 13 i=1,nvar
yp(i,kount)=y(i)
13 continue
xsav=x
endif
endif
endif
if((x+h-x2)*(x+h-x1).gt.0.) h=x2-x
call rkqs(y,dydx,nvar,x,h,eps,yscal,hdid,hnext,derivs)
if(hdid.eq.h)then
nok=nok+1
else
nbad=nbad+1
endif
if((x-x2)*(x2-x1).ge.0.)then
do 14 i=1,nvar
ystart(i)=y(i)
14 continue
if(kmax.ne.0)then
kount=kount+1
xp(kount)=x
do 15 i=1,nvar
yp(i,kount)=y(i)
15 continue
endif
return
endif
if(abs(hnext).lt.hmin) pause
*'stepsize smaller than minimum in odeint'
h=hnext
16 continue
pause 'too many steps in odeint'
return
END

SUBROUTINE bsstep(y,dydx,nv,x,htry,eps,yscal,hdid,hnext,derivs)
INTEGER nv,NMAX,KMAXX,IMAX
REAL*8 eps,hdid,hnext,htry,x,dydx(nv),y(nv),yscal(nv),SAFE1,SAFE2,
*REDMAX,REDMIN,TINY,SCALMX
PARAMETER (NMAX=50,KMAXX=8,IMAX=KMAXX+1,SAFE1=.25D0,SAFE2=.7D0,
*REDMAX=1.d-5,REDMIN=.7D0,TINY=1.d-30,SCALMX=.1D0)
CU USES derivs,mmid,pzextr
INTEGER i,iq,k,kk,km,kmax,kopt,nseq(IMAX)
REAL*8 eps1,epsold,errmax,fact,h,red,scale,work,wrkmin,xest,xnew,
*a(IMAX),alf(KMAXX,KMAXX),err(KMAXX),yerr(NMAX),ysav(NMAX),
*yseq(NMAX)
LOGICAL first,reduct
SAVE a,alf,epsold,first,kmax,kopt,nseq,xnew
EXTERNAL derivs
DATA first/.true./,epsold/-1.D0/
DATA nseq /2,4,6,8,10,12,14,16,18/
if(eps.ne.epsold)then
hnext=-1.d29
xnew=-1.d29
eps1=SAFE1*eps
a(1)=nseq(1)+1
do 11 k=1,KMAXX
a(k+1)=a(k)+nseq(k+1)
11 continue
do 13 iq=2,KMAXX
do 12 k=1,iq-1
alf(k,iq)=eps1**((a(k+1)-a(iq+1))/((a(iq+1)-a(1)+1.)*(2*k+
*1)))
12 continue
13 continue
epsold=eps
do 14 kopt=2,KMAXX-1
if(a(kopt+1).gt.a(kopt)*alf(kopt-1,kopt))goto 1
14 continue
1 kmax=kopt
endif
h=htry
do 15 i=1,nv
ysav(i)=y(i)
15 continue
if(h.ne.hnext.or.x.ne.xnew)then
first=.true.
kopt=kmax
endif
reduct=.false.
2 do 17 k=1,kmax
xnew=x+h
if(xnew.eq.x)pause 'step size underflow in bsstep'
call mmid(ysav,dydx,nv,x,h,nseq(k),yseq,derivs)
xest=(h/nseq(k))**2
call pzextr(k,xest,yseq,y,yerr,nv)
if(k.ne.1)then
errmax=TINY
do 16 i=1,nv
errmax=max(errmax,abs(yerr(i)/yscal(i)))
16 continue
errmax=errmax/eps
km=k-1
err(km)=(errmax/SAFE1)**(1./(2*km+1))
endif
if(k.ne.1.and.(k.ge.kopt-1.or.first))then
if(errmax.lt.1.)goto 4
if(k.eq.kmax.or.k.eq.kopt+1)then
red=SAFE2/err(km)
goto 3
else if(k.eq.kopt)then
if(alf(kopt-1,kopt).lt.err(km))then
red=1./err(km)
goto 3
endif
else if(kopt.eq.kmax)then
if(alf(km,kmax-1).lt.err(km))then
red=alf(km,kmax-1)*SAFE2/err(km)
goto 3
endif
else if(alf(km,kopt).lt.err(km))then
red=alf(km,kopt-1)/err(km)
goto 3
endif
endif
17 continue
3 red=min(red,REDMIN)
red=max(red,REDMAX)
h=h*red
reduct=.true.
goto 2
4 x=xnew
hdid=h
first=.false.
wrkmin=1.d35
do 18 kk=1,km
fact=max(err(kk),SCALMX)
work=fact*a(kk+1)
if(work.lt.wrkmin)then
scale=fact
wrkmin=work
kopt=kk+1
endif
18 continue
hnext=h/scale
if(kopt.ge.k.and.kopt.ne.kmax.and..not.reduct)then
fact=max(scale/alf(kopt-1,kopt),SCALMX)
if(a(kopt+1)*fact.le.wrkmin)then
hnext=h/fact
kopt=kopt+1
endif
endif
return
END

SUBROUTINE mmid(y,dydx,nvar,xs,htot,nstep,yout,derivs)
INTEGER nstep,nvar,NMAX
REAL*8 htot,xs,dydx(nvar),y(nvar),yout(nvar)
EXTERNAL derivs
PARAMETER (NMAX=50)
INTEGER i,n
REAL*8 h,h2,swap,x,ym(NMAX),yn(NMAX)
h=htot/nstep
do 11 i=1,nvar
ym(i)=y(i)
yn(i)=y(i)+h*dydx(i)
11 continue
x=xs+h
call derivs(x,yn,yout)
h2=2.*h
do 13 n=2,nstep
do 12 i=1,nvar
swap=ym(i)+h2*yout(i)
ym(i)=yn(i)
yn(i)=swap
12 continue
x=x+h
call derivs(x,yn,yout)
13 continue
do 14 i=1,nvar
yout(i)=0.5*(ym(i)+yn(i)+h*yout(i))
14 continue
return
END

SUBROUTINE pzextr(iest,xest,yest,yz,dy,nv)
INTEGER iest,nv,IMAX,NMAX
REAL*8 xest,dy(nv),yest(nv),yz(nv)
PARAMETER (IMAX=13,NMAX=50)
INTEGER j,k1
REAL*8 delta,f1,f2,q,d(NMAX),qcol(NMAX,IMAX),x(IMAX)
SAVE qcol,x
x(iest)=xest
do 11 j=1,nv
dy(j)=yest(j)
yz(j)=yest(j)
11 continue
if(iest.eq.1) then
do 12 j=1,nv
qcol(j,1)=yest(j)
12 continue
else
do 13 j=1,nv
d(j)=yest(j)
13 continue
do 15 k1=1,iest-1
delta=1./(x(iest-k1)-xest)
f1=xest*delta
f2=x(iest-k1)*delta
do 14 j=1,nv
q=qcol(j,k1)
qcol(j,k1)=dy(j)
delta=d(j)-q
dy(j)=f1*delta
d(j)=f2*delta
yz(j)=yz(j)+dy(j)
14 continue
15 continue
do 16 j=1,nv
qcol(j,iest)=dy(j)
16 continue
endif
return
END

SUBROUTINE drv(X,Y,DYDX)
IMPLICIT NONE
INTEGER NMAX
PARAMETER (NMAX=50)
REAL*8 X,DYDX(NMAX),Y(NMAX)
REAL*8 mu,nu,sigma,sigmap
REAL*8 dmudx,dnudx,dsigmadx,dsigmapdx
REAL*8 Omega,Lambda
COMMON/para/ Omega, Lambda
mu=Y(1)
nu=Y(2)
sigma=Y(3)
sigmap=Y(4)
IF(DABS(X).LT.1.0d-8) THEN
dmudx=(Omega**2*dexp(-nu)*sigma**2+sigma**2+Lambda/2.0*sigma**4)
1*x*dexp(mu)
dmudx=(Omega**2*dexp(-nu)*sigma**2-sigma**2-Lambda/2.0*sigma**4)
1*x*dexp(mu)
dsigmadx=sigmap
dsigmapdx=-dexp(mu-nu)*Omega**2*sigma+dexp(mu)*(1+Lambda*sigma**
12)*sigma
ELSE
dmudx=(Omega**2*dexp(-nu)*sigma**2+dexp(-mu)*sigmap**2+sigma**2
1+Lambda/2*sigma**4)*x*dexp(mu)-(dexp(mu)-1)/x
dnudx=(Omega**2*dexp(-nu)*sigma**2+dexp(-mu)*sigmap**2-sigma**2
1-Lambda/2*sigma**4)*x*dexp(mu)+(dexp(mu)-1)/x
dsigmadx=sigmap
dsigmapdx=((sigma**2+Lambda/2*sigma**4)*x*dexp(mu)-(dexp(mu)-1)
1/x-2/x)*sigmap-dexp(mu-nu)*Omega**2*sigma+dexp(mu)*(1+Lambda*
2sigma**2)*sigma;
ENDIF
DYDX(1)=dmudx
DYDX(2)=dnudx
DYDX(3)=dsigmadx
DYDX(4)=dsigmapdx

RETURN
END