! Here are useful F programs for Astr137   v1.3
! To use this module, compile it once (g95 -c a137.f95 ) and
! thereafter copy a137.o & a137.mod into the folder (if different)
! where you will be writing your programs. Insert   USE a137   
! after the program statement in your application. Build as 
! usual.
!
! Contents (not including private functions):
!  IMage  	Draw an image in color or greyscale in GrWin
!  IMhist  	Draw binned or cumulative histograms
!  moment  	Data moments (mean, variance, etc)
!  qsort        sort real array into ascending sequence
!  gaussdev  	Random numbers from Gaussian distribution
!  poissondev	Random number from Poisson distribution
! (Use  Fortran 95 intrinsic function  
!	call random_number(array)   for array of random #s
!				    uniform between 0 & 1 )
!  realft  	1d Fourier Transform
!  rlft2   	2d Fourier Transform
!  convlv  	convolution of two 1d functions
!  spctrm  	power spectrum of 1d data stream
!  wfits	write FITS file (uses cfitsio)

module A137	! useful programs for Astr 137
	integer, parameter, private :: dp = selected_real_kind(precision(0.0)+1)
	real(kind=dp), public, parameter :: PI = 3.14159265358979323846264338_dp
	integer, public, parameter :: erase_pen = 0, black_pen = 1, red_pen = 2, green_pen = 3, &
		blue_pen = 4, cyan_pen = 5, magenta_pen = 6, yellow_pen = 7
	interface assert_eq
		module procedure assert_eq2, assert_eq3
	end interface
	interface assert
		module procedure assert_1, assert_m
	end interface
	interface swap
		module procedure swap_r, swap_cmplx
	end interface
	interface arth
		module procedure arth_r, arth_i
	end interface
	public :: DSgraph, IMage, spectrm, convlv, four1, realft, four2, qsort, gaussdev, &
		poissondev, rlft2, moment, wfits
	private :: swap, assert_eq, assert_eq2, assert_eq3, &
		swapm, array_copy, gammln, fourrow, &
		swap_r, swap_cmplx, assert, arth, arth_r, arth_i, &
		assert_1, assert_m, zroots_unity, window
contains

recursive function zroots_unity(n,nn) result(zroots_unit)
integer, intent(in) :: n,nn
complex, dimension(nn) :: zroots_unit
integer :: k
real :: theta
!=============================================================
zroots_unit(1) = 1.0
theta = 2*PI/n
k = 1
do
	if ( k >= nn ) exit
	zroots_unit(k+1) = cmplx(cos(k*theta),sin(k*theta))
	zroots_unit(k+2:min(2*k,nn)) = zroots_unit(k+1)* &
		zroots_unit(2:min(k,nn-k))
	k = 2*k
end do
end function zroots_unity

subroutine array_copy(src,dest,n_copied,n_not_copied)
real, dimension(:), intent(in) :: src
real, dimension(:), intent(out) :: dest
integer, intent(out) :: n_copied, n_not_copied
!=============================================================
n_copied = min(size(src),size(dest))
n_not_copied = size(src)-n_copied
dest(1:n_copied) = src(1:n_copied)
end subroutine array_copy

function assert_eq2(n1,n2,string) result(assert)
character(len=*), intent(in) :: string
integer, intent(in) :: n1,n2
integer :: assert
if ( n1 == n2 ) then
	assert = n1
else
	print*,"nerror: an assert_eq fialed with this tag:",string
end if
end function assert_eq2

function assert_eq3(n1,n2,n3,string) result(assert)
character(len=*), intent(in) :: string
integer, intent(in) :: n1,n2,n3
integer :: assert
if ( n1 == n2 .and. n2 == n3 ) then
	assert = n1
else
	print*,"nerror: an assert_eq failed with this tag:",string
end if
end function assert_eq3

subroutine assert_m(n1,string)
character(len=*), intent(in) :: string
logical, intent(in), dimension(:) :: n1
if ( .not. all(n1) ) then
	print*,"nerror: an assertion failed with this tag:",string
	stop
end if
end subroutine assert_m

subroutine assert_1(n1,string)
character(len=*), intent(in) :: string
logical, intent(in) :: n1
if ( .not. n1 ) then
	print*,"nerror: an assertion failed with this tag:",string
	stop
end if
end subroutine assert_1

function arth_r(first,increment,n) result(arthval)
real(kind=dp), intent(in) :: first,increment
integer, intent(in) :: n
real(kind=dp), dimension(n) :: arthval
integer :: k,k2
real(kind=dp) :: temp
integer, parameter :: NPAR_ARTH=16, NPAR2_ARTH=8
!===============================================================
if ( n > 0 ) arthval(1) = first
if ( n <= NPAR_ARTH ) then
	do k=2,n
		arthval(k) = arthval(k-1)+increment
	end do
else
	do k=2,NPAR2_ARTH
		arthval(k) = arthval(k-1)+increment
	end do
	temp = increment*NPAR2_ARTH
	k = NPAR2_ARTH
	do
		if ( k >= n ) exit
		k2 = k+k
		arthval(k+1:min(k2,n)) = temp+arthval(1:min(k,n-k))
		temp = temp+temp
		k = k2
	end do
end if
end function arth_r

function arth_i(first,increment,n) result(arthval)
integer, intent(in) :: first,increment,n
integer, dimension(n) :: arthval
integer :: k,k2,temp
integer, parameter :: NPAR_ARTH=16, NPAR2_ARTH=8
if ( n > 0 ) arthval(1) = first
if ( n <= NPAR_ARTH ) then
	do k=2,n
		arthval(k) = arthval(k-1)+increment
	end do
else
	do k=2,NPAR2_ARTH
		arthval(k) = arthval(k-1)+increment
	end do
	temp = increment*NPAR2_ARTH
	k = NPAR2_ARTH
	do
		if ( k >= n ) exit
		k2 = k+k
		arthval(k+1:min(k2,n)) = temp+arthval(1:min(k,n-k))
		temp = temp+temp
		k = k2
	end do
end if
end function arth_i

subroutine realft(data,isign,zdata)
real, dimension(:), intent(inout) :: data
integer, intent(in) :: isign
complex, dimension(:), intent(out), optional, target :: zdata
integer :: n,ndum,nh,nq
complex, dimension(size(data)/4) :: w
complex, dimension(size(data)/4-1) :: h1,h2
complex, dimension(:), pointer :: cdata
complex :: z
real :: c1=0.5,c2
!==============================================================
n = size(data)
call assert(iand(n,n-1)==0,"n must be a power of 2 in realft")
nh = n/2
nq = n/4
if ( present(zdata) ) then
	ndum = assert_eq(n/2,size(zdata),"realft")
	cdata=>zdata
	if ( isign == 1 ) cdata = cmplx(data(1:n-1:2),data(2:n:2))
else
	allocate(cdata(n/2))
	cdata = cmplx(data(1:n-1:2),data(2:n:2))
end if
if ( isign == 1 ) then
	c2 = -0.5
	call four1(cdata,+1)
else
	c2 = 0.5
end if
w = zroots_unity(sign(n,isign),n/4)
w = cmplx(-aimag(w),real(w))
h1 = c1*(cdata(2:nq)+conjg(cdata(nh:nq+2:-1)))
h2 = c2*(cdata(2:nq)-conjg(cdata(nh:nq+2:-1)))
cdata(2:nq) = h1+w(2:nq)*h2
cdata(nh:nq+2:-1) = conjg(h1-w(2:nq)*h2)
z = cdata(1)
if ( isign == 1 ) then
	cdata(1) = cmplx(real(z)+aimag(z),real(z)-aimag(z))
else
	cdata(1) = cmplx(c1*(real(z)+aimag(z)),c1*(real(z)-aimag(z)))
	call four1(cdata,-1)
end if
if ( present(zdata) ) then
	if ( isign /= 1 ) then
		data(1:n-1:2) = real(cdata)
		data(2:n:2) = aimag(cdata)
	end if
else
	data(1:n-1:2) = real(cdata)
	data(2:n:2) = aimag(cdata)
	deallocate(cdata)
end if
end subroutine realft

subroutine moment(data,ave,adev,sdev,var,skew,curt)
real, intent(out) :: ave,adev,sdev,var,skew,curt
real, dimension(:), intent(in) :: data
integer :: n
real :: ep
real, dimension(size(data)) :: p,s
!=============================================================
n = size(data)
ave = sum(data)/n
s = data - ave
ep = sum(s)
adev = sum(abs(s))/n
p = s**2
var = sum(p)
p = p*s
skew = sum(p)
curt = sum(p*s)
var = ( var-ep**2/n )/(n-1)
sdev = sqrt(var)
if ( var /= 0.0 ) then
	skew = skew/(n*sdev**3)
	curt = curt/(n*var**2)-3.0
else
	print*,"moment: no skew or kurtosis when 0 variance"
	stop
end if
end subroutine moment

subroutine four1(data,isign)
complex, dimension(:), intent(inout) :: data
integer, intent(in) :: isign
complex, dimension(:,:), allocatable :: dat,temp
complex(kind=dp), dimension(:), allocatable :: w,wp
real(kind=dp), dimension(:), allocatable :: theta
integer :: n,m1,m2,j
!============================================================
n = size(data)
call assert(iand(n,n-1)==0,"n must be power of 2 in four1")
m1 = 2**ceiling(0.5*log(real(n))/0.693147)
m2 = n/m1
allocate(dat(m1,m2),theta(m1),w(m1),wp(m1),temp(m2,m1))
dat = reshape(data,shape(dat))
call fourrow(dat,isign)
theta = arth(0,isign,m1)*2*PI/n
wp = cmplx(-2.0_dp*sin(0.5_dp*theta)**2,sin(theta),kind=dp)
w = cmplx(1.0_dp,0.0_dp,kind=dp)
do j=2,m2
	w = w*wp+w
	dat(:,j) = dat(:,j)*w
end do
temp = transpose(dat)
call fourrow(temp,isign)
data = reshape(temp,shape(data))
deallocate(dat,w,wp,theta,temp)
end subroutine four1

subroutine rlft2(data,spec,speq,isign)
real, dimension(:,:), intent(inout) :: data
complex, dimension(:,:), intent(out) :: spec
complex, dimension(:), intent(out) :: speq
integer, intent(in) :: isign
integer :: i1,j1,nn1,nn2
real(kind=dp) :: theta
complex :: c1=(0.5,0.0),c2,h1,h2,w
complex, dimension(size(data,2)-1) :: h1a,h2a
complex(kind=dp) :: ww,wp
!============================================================
nn1 = assert_eq(size(data,1),2*size(spec,1),"rlft2: nn1")
nn2 = assert_eq(size(data,2),size(spec,2),size(speq),"rlft2: nn2")
call assert(iand((/nn1,nn2/),(/nn1,nn2/)-1)==0, &
	"dimensions must be powers of 2 in rlft2")
c2 = cmplx(0.0,-0.5*isign)
theta = 2*PI/(isign*nn1)
wp = cmplx(-2.0_dp*sin(0.5_dp*theta)**2,sin(theta))
if ( isign == 1 ) then
	spec(:,:) = cmplx(data(1:nn1:2,:),data(2:nn1:2,:))
	call four2(spec,isign)
	speq = spec(1,:)
end if
h1 = c1*(spec(1,1)+conjg(speq(1)))
h1a = c1*(spec(1,2:nn2)+conjg(speq(nn2:2:-1)))
h2 = c2*(spec(1,1)-conjg(speq(1)))
h2a = c2*(spec(1,2:nn2)-conjg(speq(nn2:2:-1)))
spec(1,1) = h1+h2
spec(1,2:nn2) = h1a+h2a
speq(1) = conjg(h1-h2)
speq(nn2:2:-1) = conjg(h1a-h2a)
ww = cmplx(1.0_dp,0.0_dp,kind=dp)
do i1=2,nn1/4+1
	j1 = nn1/2-i1+2
	ww = ww*wp+ww
	w = ww
	h1 = c1*(spec(i1,1)+conjg(spec(j1,1)))
	h1a = c1*(spec(i1,2:nn2)+conjg(spec(j1,nn2:2:-1)))
	h2 = c2*(spec(i1,1)-conjg(spec(j1,1)))
	h2a = c2*(spec(i1,2:nn2)-conjg(spec(j1,nn2:2:-1)))
	spec(i1,1) = h1+w*h2
	spec(i1,2:nn2) = h1a+w*h2a
	spec(j1,1) = conjg(h1-w*h2)
	spec(j1,nn2:2:-1) = conjg(h1a-w*h2a)
end do
if ( isign == -1 ) then
	call four2(spec,isign)
	data(1:nn1:2,:) = real(spec)/(nn1*nn2)
	data(2:nn1:2,:) = aimag(spec)/(nn1*nn2)
end if
end subroutine rlft2

subroutine four2(data,isign)
complex, dimension(:,:), intent(inout) :: data
integer, intent(in) :: isign
complex, dimension(size(data,2),size(data,1)) :: temp
!==========================================================
call fourrow(data,isign)
temp = transpose(data)
call fourrow(temp,isign)
data = transpose(temp)
end subroutine four2

subroutine convlv(data,respns,isign,conv)
real, dimension(:), intent(inout) :: data
real, dimension(:), intent(in) :: respns
integer, intent(in) :: isign
real, dimension(:), intent(out) :: conv
integer :: no2,n,m
complex, dimension(size(data)/2) :: tmpd,tmpr
!===========================================================
n = size(data)
m = size(respns)
call assert(iand(n,n-1)==0,"n must be a power of 2 in convlv")
call assert(modulo(m,2)== 1,"m must be odd in convlv")
conv(1:m) = respns(:)
conv(n-(m-3)/2:n) = conv((m+3)/2:m)
conv((m+3)/2:n-(m-1)/2) = 0.0
no2 = n/2
call realft(data,1,tmpd)
call realft(conv,1,tmpr)
if ( isign == 1 ) then
	tmpr(1) = cmplx(real(tmpd(1))*real(tmpr(1))/no2, &
		aimag(tmpd(1))*aimag(tmpr(1))/no2)
	tmpr(2:) = tmpd(2:)*tmpr(2:)/no2
else if ( isign == -1 ) then
	if ( any(abs(tmpr(2:)) == 0.0 ) .or. real(tmpr(1)) == 0.0 &
		.or. aimag(tmpr(1)) == 0.0 ) print*,"deconvoling at response 0 in convlv"
	tmpr(1) = cmplx(real(tmpd(1))/real(tmpr(1))/no2, &
		aimag(tmpr(1))/aimag(tmpr(1))/no2)
	tmpr(2:) = tmpd(2:)/tmpr(2:)/no2
else
	print*,"no meaning for isign in convlv"
end if
call realft(conv,-1,tmpr)
end subroutine convlv

subroutine spectrm(p,k,ovrlap,unit,n_window)
real, dimension(:), intent(inout) :: p
integer, intent(in) :: k
logical, intent(in) :: ovrlap
integer, optional, intent(in) :: n_window,unit
integer :: j,joff,joffn,kk,m,m4,m43,m44,mm,iunit,nn_window
real :: den,facm,facp,sumw
real, dimension(2*size(p)) :: w
real, dimension(4*size(p)) :: w1
real, dimension(size(p)) :: w2
complex, dimension(2*size(p)) :: cw1
!========================================================================================
m = size(p)
if ( present(n_window) ) then
	nn_window = n_window
else
	nn_window = 1
end if
if ( present(unit) ) then
	iunit = unit
else
	iunit = 9
end if
mm = m+m
m4 = mm+mm
m44 = m4+4
m43 = m4+3
den = 0.0
facm = m
facp = 1.0/m
w1(1:mm) = window(arth(1,1,mm),facm,facp,nn_window)
sumw = dot_product(w1(1:mm),w1(1:mm))
p = 0.0
if ( ovrlap ) read(unit=iunit,fmt=*) (w2(j),j=1,m)
do kk=1,k
	do joff=-1,0,1
		if ( ovrlap ) then
			w1(joff+2:joff+mm:2) = w2(1:m)
			read(unit=iunit,fmt=*) (w2(j),j=1,m)
			joffn = joff+mm
			w1(joff+2:joffn+mm:2) = w2(1:m)
		else
			read(unit=iunit,fmt=*) (w1(j),j=joff+2,m4,2)
		end if
	end do
	w = window(arth(1,1,mm),facm,facp,nn_window)
	w1(2:m4:2) = w1(2:m4:2)*w
	w1(1:m4:2) = w1(1:m4:2)*w
	cw1(1:mm) = cmplx(w1(1:m4:2),w1(2:m4:2))
	call four1(cw1(1:mm),1)
	w1(1:m4:2) = real(cw1(1:mm))
	w1(2:m4:2) = aimag(cw1(1:mm))
	p(1) = p(1)+w(1)**2+w1(2)**2
	p(2:m) = p(2:m)+w1(4:2*m:2)**2+w1(3:2*m-1:2)**2 + &
		w1(m44-4:m44-2*m:-2)**2+w1(m43-4:m43-2*m:-2)**2
	den = den+sumw
end do
p = p/(m4*den)
end subroutine spectrm

function window(j,facm,facp,nn_window) result(wind)
integer, dimension(:), intent(in) :: j
integer, intent(in) :: nn_window
real, intent(in) :: facm,facp
real, dimension(size(j)) :: wind
select case(nn_window)
	case(1)
		wind(j) = 1.0-abs((j-1)-facm)*facp
	case(2)
		wind(j) = 1.0
	case(3)
		wind(j) = 1.0-(((j-1)-facm)*facp)**2
	case default
		print*,"unimplemented window function in spctrm"
end select
end function window

subroutine fourrow(data,isign)
complex, dimension(:,:), intent(inout) :: data
integer, intent(in) :: isign
integer :: n,i,istep,j,m,mmax,n2
real(kind=dp) :: theta
complex, dimension(size(data,1)) :: temp
complex(kind=dp) :: w,wp
complex :: ws
!===========================================================
n = size(data,2)
n2 = n/2
j = n2
do i=1,n-2
	if ( j > i ) call swap(data(:,j+1),data(:,i+1))
	m = n2
	do 
		if ( m < 2 .or. j < m ) exit
		j = j-m
		m = m/2
	end do
	j = j+m
end do
mmax = 1
do
	if ( n <= mmax ) exit
	istep = 2*mmax
	theta = PI/(isign*mmax)
	wp = cmplx(-2.0_dp*sin(0.5_dp*theta)**2,sin(theta),kind=dp)
	w = cmplx(1.0_dp,0.0_dp,kind=dp)
	do m=1,mmax
		ws = w
		do i=m,n,istep
			j = i+mmax
			temp = ws*data(:,j)
			data(:,j) = data(:,i)-temp
			data(:,i) = data(:,i)+temp
		end do
		w = w*wp+w
	end do
	mmax = istep
end do
end subroutine fourrow

subroutine gaussdev(harvest)
! return array of random numbers distributed in Gaussian
! of mean 0 and sigma 1
real, dimension(:), intent(out) :: harvest
real, allocatable, dimension(:), save :: g
real, dimension(size(harvest)) :: rsq,v1,v2
integer :: n,ng,nn,m
integer, save :: last_allocated = 0
logical, save :: gaus_stored = .false.
logical, dimension(size(harvest)) :: mask
!==============================================================
n = size(harvest)
if ( n /= last_allocated ) then
	if ( last_allocated /= 0 ) deallocate(g)
	allocate(g(n))
	last_allocated = n
	gaus_stored = .false.
end if
if ( gaus_stored ) then
	harvest = g
	gaus_stored = .false.
else
	ng = 1
	do
		if ( ng > n ) exit
		call random_number(v1(ng:n))
		call random_number(v2(ng:n))
		v1(ng:n) = 2.0*v1(ng:n)-1.0
		v2(ng:n) = 2.0*v2(ng:n)-1.0
		rsq(ng:n) = v1(ng:n)**2 + v2(ng:n)**2
		mask(ng:n) = ( rsq(ng:n) > 0.0 .and. rsq(ng:n) < 1.0 )
		call array_copy(pack(v1(ng:n),mask(ng:n)),v1(ng:),nn,m)
		v2(ng:ng+nn-1) = pack(v2(ng:n),mask(ng:n))
		rsq(ng:ng+nn-1) = pack(rsq(ng:n),mask(ng:n))
		ng = ng+nn
	end do
	rsq = sqrt(-2.0*log(rsq)/rsq)
	harvest = v1*rsq
	g = v2*rsq
	gaus_stored = .true.
end if
end subroutine gaussdev

function gammln(xx) result(gamln)
real, intent(in) :: xx
real :: gamln
real(kind=dp) :: tmp,x
real(kind=dp), parameter :: stp = 2.5066282746310005_dp
real(kind=dp), dimension(6), parameter :: coef = (/ 76.18009172947146_dp, &
	-86.50532032941677_dp, 24.01409824083091_dp, -1.231739572450155_dp, &
	0.1208650973866179e-2_dp, -0.5395239384953e-5_dp /)
!========================================================
x = xx
tmp = x+5.5_dp
tmp = (x+0.5_dp)*log(tmp)-tmp
gamln = tmp+log(stp*(1.000000000190015_dp + &
	sum(coef(:)/arth(x+1.0_dp,1.0_dp,size(coef))))/x)
end function gammln

subroutine poissondev(xm,poidev)
! return one Poisson deviate. Note: F cannot do this in function
! because it assumes that functions are PURE, i.e. cannot SAVE anything
real, intent(in) :: xm
real, intent(out) :: poidev
real :: em,harvest,t,y
real, save :: alxm, g, sq, oldm=-1.0
!=======================================================
if ( xm < 12.0 ) then
	if ( xm /= oldm ) then
		oldm = xm
		g = exp(-xm)
	end if
	em = -1
	t = 1.0
	do
		em = em+1.0
		call random_number(harvest)
		t = t*harvest
		if ( t <= g ) exit
	end do
else
	if ( xm /= oldm ) then
		oldm = xm
		sq = sqrt(2.0*xm)
		alxm = log(xm)
		g = xm*alxm-gammln(xm+1.0)
	end if
	do
		do
			call random_number(harvest)
			y = tan(PI*harvest)
			em = sq*y+xm
			if ( em >= 0.0 ) exit
		end do
		em = int(em)
		t = 0.9*(1.0+y**2)*exp(em*alxm-gammln(em+1.0)-g)
		call random_number(harvest)
		if ( harvest <= t ) exit
	end do
end if
poidev = em
end subroutine poissondev

subroutine swap_r(a,b)
real, intent(inout) :: a,b
real :: dum
dum = a
a = b
b = dum
end subroutine swap_r

subroutine swap_cmplx(a,b)
complex, dimension(:), intent(inout) :: a,b
complex, dimension(size(a,1)) :: dum
dum = a
a = b
b = dum
end subroutine swap_cmplx

subroutine swapm(a,b,mask)
real, intent(inout) :: a,b
logical, intent(in) :: mask
real :: swp
if ( mask ) then
	swp = a
	a = b
	b = swp
end if
end subroutine swapm

subroutine qsort(arr)	
! sort real array into ascending order, overwriting the input
real, dimension(:), intent(inout) :: arr
integer, parameter :: NN=15, NSTACK=50
real :: a
integer :: n,k,i,j,jstack,l,r
integer, dimension(NSTACK) :: istack
n = size(arr)
jstack = 0
l = 1
r = n
do
	if ( r-l < NN ) then
		do j=l+1,r
			a = arr(j)
			do i=j-1,l,-1
				if ( arr(i) <= a ) exit
				arr(i+1) = arr(i)
			end do
			arr(i+1) = a
		end do
		if ( jstack == 0 ) RETURN
		r = istack(jstack)
		l = istack(jstack-1)
		jstack = jstack-2
	else
		k = (l+r)/2
		call swap(arr(k),arr(l+1))
		call swapm(arr(l),arr(r),arr(l)>arr(r))
		call swapm(arr(l+1),arr(r),arr(l+1)>arr(r))
		call swapm(arr(l),arr(l+1),arr(l)>arr(l+1))
		i = l+1
		j = r
		a = arr(l+1)
		do
			do
				i = i+1
				if ( arr(i) >= a ) exit
			end do
			do
				j = j-1
				if ( arr(j) <= a ) exit
			end do
			if ( j < i ) exit
			call swap(arr(i),arr(j))
		end do
		arr(l+1) = arr(j)
		arr(j) = a
		jstack = jstack+2
		if ( jstack > NSTACK ) then
			print*,"sort: NSTACK too small"
			stop
		end if
		if ( r-i+1 >= j-l ) then
			istack(jstack) = r
			istack(jstack-1) = i
			r = j-1
		else
			istack(jstack) = j-1
			istack(jstack-1) = l
			l = i
		end if
	end if
end do
end subroutine qsort

subroutine IMage(type,im,cap)
real, dimension(:,:), intent(in) :: im
character(len=*), dimension(3), intent(in), optional :: cap
character(len=*), intent(in) :: type
!=================================================================================
if ( present(cap) ) then
	if ( len_trim(cap(1)) /= 0 ) call name(cap(1),"x")
	if ( len_trim(cap(2)) /= 0 ) call name(cap(2),"y")
	if ( len_trim(cap(3)) /= 0 ) call titlin(cap(3),2)
end if
call setvlt(type)
call qplclr(im,size(im,1),size(im,2))
end subroutine IMage

subroutine DSgraph(device)
character(len=*), intent(in) :: device
!=================================================================================
call setpag("USAL")
if ( device == "xwin" .or. device(1:3) /= "eps" ) then
	call metafl(device)
else if ( index(device,"P") == 0 ) then
	call metafl(device(1:3))
else 
	call setpag("USAP")
	call metafl(device(1:3))
end if
call scrmod("reverse")
call disini()
call axslen(1500,1500)
end subroutine DSgraph

! subroutine IMshow(type,f,i1in,i2in,i3in,i4in,fmin,fmax,t)
! ! Display image on graphics device
! ! type is colortable: 1= greyscale, 2= rainbow, 3 = heat, etc
! ! 	(if -ve then add to existing page)
! ! f: array dimensioned containing values to plot
! ! 4 boundaries: optional
! ! 	i1: bottom left x, i2: bottom left y, i3: top right x, i4: top right y
! ! fmin: minimum intensity to plot, values below are ignored (optional)
! ! fmax: maximum intensity, values above are ignored (optional)
! ! t: rotation/shear coefficients (optional
! integer, intent(in) :: type
! real, dimension(:,:), intent(in) :: f
! real, dimension(6), intent(in), optional :: t
! real, intent(in), optional :: i1in,i2in,i3in,i4in,fmin,fmax
! real :: i1,i2,i3,i4,fmn,fmx
! real, parameter :: contra = 1.0, bright = 0.0
! real, dimension(2) :: gl, gr, gg, gb
! real, dimension(9) :: rl, rr, rg, rb
! real, dimension(5) :: hl, hr, hg, hb
! real, dimension(10) :: wl, wr, wg, wb
! real, dimension(20) :: al, ar, ag, ab
! real, dimension(6) :: tr,ts = (/0, 1, 0, 0, 0, 1/)
! real :: d, vpx1, vpx2, vpy1, vpy2
! integer :: lut,mxi,mxj
! !========================================================
! mxi = size(f,1)
! mxj = size(f,2)
! if ( present(t) ) then
! 	tr = t
! else
! 	tr = ts
! end if
! if ( present(i1in) ) then
! 	i1 = i1in
! else
! 	i1 = 1.0
! end if
! if ( present(i2in) ) then
! 	i2 = i2in
! else
! 	i2 = size(f,1)
! end if
! if ( present(i3in) ) then
! 	i3 = i3in
! else
! 	i3 = 1.0
! end if
! if ( present(i4in) ) then
! 	i4 = i4in
! else
! 	i4 = size(f,2)
! end if
! if ( present(fmin) ) then
! 	fmn = fmin
! else
! 	fmn = minval(f)
! end if
! if ( present(fmax) ) then
! 	fmx = fmax
! else
! 	fmx = maxval(f)
! end if
! if ( type > 0 ) then
! 	call pgpage()
! 	call pgsvp(0.0, 1.0, 0.0, 1.0)
! 	lut = type
! else
! 	lut = -type
! end if
! call pgqvp(1, vpx1, vpx2, vpy1, vpy2)
! d = min(vpx2-vpx1, vpy2-vpy1)/40.0
! vpx1 = vpx1 + 5.0*d
! vpx2 = vpx2 - 2.0*d
! vpy1 = vpy1 + 8.0*d -1.2
! vpy2 = vpy2 - 2.0*d
! call pgvsiz(vpx1, vpx2, vpy1, vpy2)
! gl = (/0.0, 1.0/)
! gr = gl
! gg = gl
! gb = gl
! rl = (/-0.5, 0.0, 0.17, 0.33, 0.50, 0.67, 0.83, 1.0, 1.7/)
! rr = (/ 0.0, 0.0,  0.0,  0.0,  0.6,  1.0,  1.0, 1.0, 1.0/)
! rg = (/ 0.0, 0.0,  0.0,  1.0,  1.0,  1.0,  0.6, 0.0, 1.0/)
! rb = (/ 0.0, 0.3,  0.8,  1.0,  0.3,  0.0,  0.0, 0.0, 1.0/)
! hl = (/0.0, 0.2, 0.4, 0.6, 1.0/)
! hr = (/0.0, 0.5, 1.0, 1.0, 1.0/)
! hg = (/0.0, 0.0, 0.5, 1.0, 1.0/)
! hb = (/0.0, 0.0, 0.0, 0.3, 1.0/)
! wl = (/0.0, 0.5, 0.5, 0.7, 0.7, 0.85, 0.85, 0.95, 0.95, 1.0/)
! wr = (/0.0, 1.0, 0.0, 0.0, 0.3,  0.8,  0.3,  1.0,  1.0, 1.0/)
! wg = (/0.0, 0.5, 0.4, 1.0, 0.0,  0.0,  0.2,  0.7,  1.0, 1.0/)
! wb = (/0.0, 0.0, 0.0, 0.0, 0.4,  1.0,  0.0,  0.0, 0.95, 1.0/)
! al = (/0.0, 0.1, 0.1, 0.2, 0.2, 0.3, 0.3, 0.4, 0.4, 0.5, &
!       0.5, 0.6, 0.6, 0.7, 0.7, 0.8, 0.8, 0.9, 0.9, 1.0/)
! ar = (/0.0, 0.0, 0.3, 0.3, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0, &
!       0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0/)
! ag = (/0.0, 0.0, 0.3, 0.3, 0.0, 0.0, 0.0, 0.0, 0.8, 0.8, &
!       0.6, 0.6, 1.0, 1.0, 1.0, 1.0, 0.8, 0.8, 0.0, 0.0/)
! ab = (/0.0, 0.0, 0.3, 0.3, 0.7, 0.7, 0.7, 0.7, 0.9, 0.9, &
!       0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)
! select case (lut)
! case (1) ! -- gray scale
! 	call pgctab(gl, gr, gg, gb, 2, contra, bright)
! case (2) ! -- rainbow
! 	call pgctab(rl, rr, rg, rb, 9, contra, bright)
! case (3) ! -- heat
! 	call pgctab(hl, hr, hg, hb, 5, contra, bright)
! case (4) ! -- iraf
! 	call pgctab(wl, wr, wg, wb, 10, contra, bright)
! case (5) ! -- aips
! 	call pgctab(al, ar, ag, ab, 20, contra, bright)
! end select
! call pgwnad(i1,i2,i3,i4)
! call pgimag(f,mxi,mxj,1,mxi,1,mxj,fmn,fmx,tr)
! if ( type < 0 ) then
! 	call pgwedg("ti", 2.5, 5.0, fmn, fmx, "pixel value")
! 	call pgbox("bctsi",0.0,0,"bctsi",0.0,0)
! else
! 	call pgwedg("ri", 2.5, 5.0, fmn, fmx, "pixel value")
! 	call pgbox("bcntsi",0.0,0,"bcntsiv",0.0,0)
! end if
! 
! end subroutine IMshow

! subroutine IMhist(f,dmin,dmax,bins) ! optional last 3 args
! ! Draw binned or cumulative histogram on graphics device
! ! f: 2d array, each row to histogram. If array is 1d (or 3d), then use
! !	call Imhist(reshape(f, /( nx*ny*nz,1 )/), ...)
! ! successive distributions in f are drawn in different colors
! ! dmin,dmax: only data between these limits are histogramed. If not
! !	provided then use entire data range
! ! bins: if this 2d real array is provided then its length gives # bins to
! !	fill, plot, & return. Max of 400 bins.  If not provided then plot a
! !	cumulative histogram (i.e. no binning)
! real, dimension(:,:), intent(in) :: f
! real, optional, intent(in) :: dmin,dmax
! real :: dmn,dmx
! real, dimension(:,:), optional, intent(inout) :: bins
! real, dimension(:), allocatable :: f1,y
! integer :: n,i,j,lo,hi,nb,nl,l
! ! call pgsci(l) = 2,... : rR, G, B, cyan, magenta, yellow
! real :: s
! !=========================================================
! if ( present(dmin) ) then
! 	dmn = dmin
! else
! 	dmn = minval(f)
! end if
! if ( present(dmax) ) then
! 	dmx = dmax
! else
! 	dmx = maxval(f)
! end if
! n = size(f,1)
! nl = size(f,2)
! if ( present(bins) ) then ! regular binned histogram
! 	nb = size(bins,1)
! 	allocate(y(1:nb))
! 	do l=1,nl
! 		bins(:,l) = 0.0
! 		do i=1,n
! 			j = nb*(f(i,l)-dmn)/(dmx-dmn)
! 			if ( j >= 1 .and. j <= min(nb,400) ) bins(j,l) = bins(j,l)+1
! 		end do
! 		if ( l == 1 ) then
! 			call pgsci(1)
! 			call pgenv(dmn,dmx,0.0,maxval(bins(1:nb,l),l),0,0)
! 			call pgslw(3)
! 			s = (dmx-dmn)/nb
! 			y(1:nb) = (/ (dmn+s*(i-1),i=1,nb) /)
! 		else
! 			call pgsci(l)
! 		end if
! 		call pgbin(nb,y(1:l),bins(:,l),.false.)
! 	end do
! 	deallocate(y)
! 	call pgslw(1)
! else ! cumulative histogram
! 	do l=1,nl
! 		allocate(f1(1:n))
! 		f1 = f(1:n,l)
! 		call qsort(f1)
! 		if ( l == 1 ) then
! 			do lo=1,n
! 				if ( f1(lo) >= dmn ) exit
! 			end do
! 			do hi=n,i,-1
! 				if ( f1(hi) <= dmx ) exit
! 			end do
! 			call pgsci(1)
! 			call pgenv(dmn,dmx,0.0,1.01,0,0)
! 		else
! 			call pgsci(l)
! 		end if
! 		call pgbin(hi-lo+1,f1(lo:hi),(/ (real(i)/n,i=lo,hi) /),.false.)
! 		deallocate(f1)
! 	end do
! end if
! call pgsci(1)
! end subroutine IMhist

subroutine wfits(filename,a)
character(len=*), intent(in) :: filename
real, dimension(:,:), intent(in) :: a
logical :: simple,existingfile
integer :: status,unit,blocksize,bitpix,fpixel,nx,ny,naxis,readwrite = 1
integer, dimension(3) :: axes
!=============================================================================
status = 0
call ftgiou(unit,status)
call ftinit(unit,filename,blocksize,status)
existingfile = ( status == 105 )
if ( existingfile ) then
	status = 0
	call ftfiou(unit,status)
	call ftgiou(unit,status)
	call ftopen(unit,filename,readwrite,blocksize,status)
end if
nx = size(a,1)
ny = size(a,2)
simple = .true.
bitpix = -32
axes(3) = 0
if ( ny > 0 ) then
	naxis = 2
	axes(2) = ny
else 
	naxis = 1
	axes(2) = 0
end if
axes(1) = nx
if ( .not.existingfile ) call ftphpr(unit,simple,bitpix,naxis,axes,0,1,.false.,status)
fpixel = 1
call ftppre(unit,1,fpixel,nx*max(1,ny),a,status)
call ftclos(unit,status)
call ftfiou(unit,status)
end subroutine wfits
end module a137