The effect of cooling on time dependent behaviour of accretion flows around black holes
Sandip K. Chakrabarti$^{1,2}$
K. Acharyya$^2$ and D. Molteni$^3$
We present the results of several numerical simulations of two dimensional axi-symmetric accretion flows around black holes using Smoothed Particle Hydrodynamics (SPH) in the presence of cooling effects. We consider both stellar black holes and super-massive black holes. We observe that due to {both radial and vertical oscillation} of shock waves in the accretion flow, the luminosity and average thermal energy content of the inner disk exhibit very interesting behaviour. When power density spectra are taken, quasi-periodic variabilities are seen at a few Hz and also occasionally at hundreds of Hz for stellar black holes. For super-massive black holes, the time scale of the oscillations ranges from hours to weeks. The power density spectra have a flat top behavior with average {rms} amplitude of a few percent and a broken power-law behavior. The break frequency is generally found to be close to the frequency where the shock oscillates.
CHAKRA_SKC.TEX
\documentstyle[referee]{l-aa}
\input psfig.tex
\def\lsim{\lower.5ex\hbox{$\; \buildrel < \over \sim \;$}}
\def\gsim{\lower.5ex\hbox{$\; \buildrel > \over \sim \;$}}
\def \simeq{\lower.3ex\hbox{$\; \buildrel \sim \over - \;$}}
\def\ch{\lower-0.55ex\hbox{--}\kern-0.55em{\lower0.15ex\hbox{$h$}}}
\def\lh{\lower-0.55ex\hbox{--}\kern-0.55em{\lower0.15ex\hbox{$\lambda$}}}
\font
\capb=cmbx8
\font
\ita=cmmi9
\font
\sla=cmss9
\font
\capt=cmti8.tfm
\renewcommand {\baselinestretch}{1.5}
\begin{document}
\title{The effect of cooling on time dependent behaviour of accretion flows around black holes}
\author{Sandip K. Chakrabarti$^{1,2}$, K. Acharyya$^2$ and D. Molteni$^3$}
\institute {$^1$ S. N. Bose National Centre for Basic Sciences, Salt Lake, Kolkata 700098, India\\
$^2$ Centre for Space Physics, Chalantika 43, Garia Station Rd., Kolkata, 700084, India\\ 
$^3$ Dipartimento di Fisica e Tecnologie Relative, Viale delle Scienze, 90128 Palermo Italy\\
e-mails: chakraba@bose.res.in, space\_phys@vsnl.com, molteni@unipa.it}
\offprints{S. K. Chakrabarti {\it chakraba@boson.bose.res.in}}
\date{Received ; accepted , }
\maketitle
\markboth{  }{}

\begin{abstract}

We present the results of several numerical simulations of two dimensional 
axi-symmetric accretion flows around black holes using Smoothed Particle Hydrodynamics (SPH)
in the presence of cooling effects. We consider both stellar black holes and 
super-massive black holes. We observe that due to {\it both radial and vertical
oscillation} of shock waves in the accretion flow, the luminosity and average thermal energy 
content of the inner disk exhibit very interesting behaviour. When power 
density spectra are taken, quasi-periodic variabilities are seen at a few 
Hz and also occasionally at hundreds of Hz for stellar black holes. 
For super-massive black holes, the time scale of the oscillations
ranges from hours to weeks. The power density spectra have a flat 
top behavior with average {\it rms} amplitude of a few percent 
and a broken power-law behavior. The break frequency is generally found to be
close to the frequency where the shock oscillates.

\end{abstract}

\section{Introduction}

X-rays emitted by galactic black hole candidates are known to exhibit quasi-periodic variations (QPVs).
Today, there are almost a dozen confirmed stellar mass black hole candidates for which 
QPVs are regularly observed, some at normal frequencies ($0.1-10$Hz) and some with
frequencies of hundreds of Hz. While it is conceptually easier to explain such variations when there is a 
hard surface, as on a neutron star (e.g., Shaham, 1987), in black hole candidates they are
difficult to understand as there are no surfaces from which a perturbation could be reflected back
so as to produce a sustained oscillation of significant amplitude. 

In a very important work, Langer, Chanmugam and Shaviv (1981) first suggested that 
in the presence of cooling effects standing shocks on a white dwarf surface can 
oscillate and in fact the oscillation may also subside when 
cooling is increased (Chanmugam, Langer, and Shaviv, 1985). While a black hole does not have a 
hard surface like neutron stars or white dwarfs, the centrifugal force in accretion matter
close to a black hole can become strong enough to actually slow down matter and form 
centrifugal-pressure-supported axisymmetric standing shocks 
(Chakrabarti, 1989) depending on specific energy and angular momentum 
(this region will be referred to as CENBOL, centrifugal pressure dominated boundary layer). 
In the particle dynamics picture, the appearance of a discontinuity can be argued this way: the 
centrifugal force $F_{cen} \propto \lambda^2/r^3$ grows more rapidly 
(for constant or almost constant specific angular momentum $\lambda$)
than the gravitational force $F_{G} \propto GM/r^2$ and they become comparable at a distance of 
$\lambda \sim (GMr)^{1/2}$, where $M$ is the mass of the central compact object.
In fluid dynamics, because of the addition of thermal pressure  and/or radiation  pressure, 
such a balance takes place farther out where a shock may develop. 
The almost constancy of angular momentum comes into being
primarily because, as has been shown by extensive work done earlier (Chakrabarti, 1996a), 
the viscosity timescale $t_{visc} \sim r^2/\nu_s$ is much longer than the
the infall time scale $r^{3/2}/(GM)^{1/2}$, $\nu_s$ being the kinematic viscosity. 
It was shown by Molteni, Sponholz \& Chakrabarti (1996, hereafter referred 
to as MSC96) and Ryu, Chakrabarti, \& Molteni (1997, hereafter 
RCM 1997) that these shocks may also oscillate in the presence of cooling. 
These oscillations are explained to be due to a resonance between the dominant 
cooling time scale and the infall time scale, and the cooling could be due to 
thermal/non-thermal radiative effects (as in MSC96), or to
dynamical cooling due to outflows (as in RCM97) or both. This shock
forms in low-angular-momentum (in the sense mentioned), quasi-spherical
flows which may or may not surround a
Keplerian disk located in the equatorial plane. In case the latter disk is present
and the soft photons emitted are intercepted, inverse-Comptonized and re-emitted by the 
post-shock  region, as in the two-component advective flow model (Chakrabarti
\& Titarchuk, 1995; hereafter CT95), then, as the shock oscillates, 
the intensity of hard X-rays is also modulated at the same frequency. 

If our shock oscillation solution is correct, then a large number of simple 
predictions could be made: (a) the harder radiation
coming from the hot, post-shock region would primarily participate in QPVs while the softer
radiation would have less power. This was verified by Chakrabarti \& Manickam 
(2000, hereafter CM00) and Rao et al. (2000). (b) the stronger a
shock is, the stronger is the thermally driven outflow rate. Thus, 
the outflow from the post-shock region depends on the shock strength, 
i.e., the spectral states (Chakrabarti, 1999). 
This has also been tested by various observations (e.g., 
Dhawan et al., 2000). (c) In case of outflows, the electron density in 
the post-shock region is reduced and the Comptonized spectrum should be softened.
Chakrabarti et al. (2003) showed that 
indeed the presence/absence of outflows does change the spectral slope 
in the hard X-ray region. These observations vindicate that the
Comptonized photon, and the major part of the outflows both
come out of the post-shock region. Thus, a thorough study of 
the behavior of CENBOL, especially when it is time dependent, is needed.

In the present paper, we carry out numerical simulation of thick, 
inviscid, rotating axisymmetric flows, in the presence of 
cooling effects for a large region of the parameter space. 
The hot electrons in the CENBOL
reprocess soft photons and generally enhance their
energy before these photons escape. The degree by which the 
energy is transferred depends on the so-called 
{\it  Compton y parameter} (Rybicki \& Lightman, 1979),
$$
y=\tau_0 (4 \Theta + 16 \Theta^2) {\rm max(1,}\tau_0)
\eqno{(1)}
$$
where $\tau_0$ is the optical depth of the CENBOL, 
and $\Theta=kT_e/m_e c^2$ is dimensionless electron temperature,
($k$, $T_e$, $m_e$ and $c$ being the Boltzmann constant, electron temperature,
electron mass and the velocity of light respectively).
If $y<<1$, nothing much happens to the photons, and for $y>>1$
the saturated Comptonization causes a Wien distribution to form 
at a high energy. The most important case is where $y$ is around $1$. In this
latter case, the energy loss from electrons to photons becomes important. To calculate this
cooling effect in a time-dependent simulation is not easy, since ideally, at 
every time step, one has to carry out a Monte-Carlo simulation of the 
configuration of the flow. Here we use a shortcut to the process. 
To include the effect of Comptonization, we borrow the
results of a Monte Carlo simulation by Dermer, Liang \& Canfield (1991, hereafter DLC91)
where an enhancement factor $\eta$ 
(the ratio of up-scattered flux to injected flux) has been calculated 
for various geometries of the medium. 
We  define an effective accretion rate ${\dot m}_c \sim {\dot m}_b / \eta$ 
such that the cooling due to inverse Comptonization per gram of matter in the post-shock region 
with ${\dot m}_c$ is the same as the bremsstrahlung cooling with accretion rate ${\dot m}_b$.
Essentially, we use bremsstrahlung cooling with an accretion rate enhanced by a factor  $\eta$
whose value depends on $\tau_0$ (DLC91; CT95). This way, our basic 
questions, namely, whether the shock exhibits an oscillation in the presence
of cooling, and whether such oscillation shows quasi-periodicity are answered.
It so happens that the power density spectra (PDS) we obtained have similar 
characteristics as of the PDSs of light curves observed 
in black hole candidates with quasi-periodic oscillations (QPOs), so we shall use the terms QPVs and QPOs interchangeably. 
As discussed in MSC96, the cause of the
QPOs or QPVs is generic, namely, rough agreement between the cooling time and the infall time. So
we believe that the exact cooling process is immaterial -- as long as the two problems have a similar cooling
time-scale $E/{\dot E}$, where $E$ is the energy density and ${\dot E}$ is the cooling rate, 
they would produce the same oscillations. Hence PDS should continue to show QPVs 
even when more rigorous treatment of Comptonization and other coolings is included. In future,
we plan to consider Comptonization more rigorously.

Unlike MSC96, where symmetry around the equatorial plane was assumed
(and thus the simulation was carried out only in one quadrant), we consider both halves 
above and below the equatorial plane. In this way, allowance is made not only for radial 
oscillations of the shocks  but also for vertical oscillations.  Earlier, in a 
simulation of a non-shock flow it was observed that the outgoing winds could interact with the 
inflow and the interface might undulate (Molteni et al. 2001). So we expect that shocks 
should also have vertical oscillations in presence of winds. In the present simulations,
we observe that the shocks do oscillate both vertically and horizontally for a
range of inflow parameters, and the PDS of the light curve obtained by
our numerical simulation does have characteristics similar to the PDS of the observed light curves
of the black hole candidates. The PDS of a numerically simulated light 
curve has a flat top up to a break frequency, and an {\it  rms} 
amplitude in the range  of a few percent depending on the accretion rate of the flow.
Considering that an axisymmetric shock fundamentally separates a flow 
into a fast moving (lower density) supersonic and a slow moving 
(higher density) sub-sonic region, it should not be surprising that its oscillation
frequency will  generally lie at or near the break frequency in the power law.
Of course, shocks are not infinitesimally thin. Neither is it true that
they exhibit only a single mode of oscillation. Thus broad and often multiple peaks
are expected. This is precisely what we see. What is more, since the inner sonic point 
also separates the flow into two parts, namely, one with the disk-like behavior and the  one with free-fall
behavior, there is often a weak peak in the PDS at some hundreds of Hz, perhaps corresponding 
to oscillation of this inner region. We shall show these as well. 

In the next section, we present the basic setup of our numerical simulation.
In Sect. 3 we present some results of the simulation for galactic and extra-galactic 
black holes and analyze them to demonstrate shock oscillations. From the PDS of the 
`simulated' light curves, we show that there are sharp peaks at certain
frequencies only. The PDS seems to have the familiar 
flat top shape with a break. In Sect. 4, we discuss possible cause of the
behaviour of the shocks. Finally, in Sect. 5 we draw our conclusions. 

\section{Setup of the numerical simulations} 

We are interested in simulating the behavior of an axisymmetric, inviscid flow
in the presence of cooling effects. Unlike an ordinary star, a black
hole is capable of accreting matter even in the absence of viscosity 
when the pressure is large enough. However, if the specific angular 
momentum is less than the marginally stable value ($3.67 GM_1/c$ for 
a Schwarzschild black hole of mass $M_1$), even a cool, non-viscid
gas can enter a black hole without any problem. We chose the 
specific angular momentum of the flow $\lambda$ to be lower compared to this value
everywhere. This makes $\lambda < (G M_1 r)^{1/2}$, the Keplerian value, 
everywhere in the flow. We call such a flow a `sub-Keplerian' flow. The motivation for
studying such a flow stems from the fact that accretion processes 
into a black hole are necessarily transonic, and the flow has to be 
sub-Keplerian near the horizon. Furthermore, a flow can have a 
substantial amount of sub-Keplerian matter itself because of
accretion of winds from the companion. The wind of velocity $v_w \sim (2GM_2/R_2)^{1/2}$, 
having `zero' angular momentum with respect to the
mass-lossing star of mass $M_2$ and radius $R_2$ has an angular momentum of 
$l_w \lsim v_w R_a$, where $R_a \sim 10^{11}$cm is the radius of the accretion disk. For typical
values: $M_2\sim 2 \times 10^{33}$gm, $R_2\sim 10^{11}$cm and we get,
in units of $2GM_1/c$, $l_w \lsim 3.2$ with $M_1\sim 10M_\odot$. Hence,
our choice of $1.75$ which is averaged over the entire 
flow (thus could be much smaller than $l_w$) is justified. 
Recent observations do indicate the presence of both Keplerian 
and sub-Keplerian components (Smith, Heindl \& Swank, 2002) in half a dozen
black hole components.

It was shown by Chakrabarti (1989, 1996a) that in the large region 
of the parameter space spanned by the specific angular momentum $\lambda$ 
and energy ${\cal E}$ a stable standing shock may form in accretion whenever 
the Rankine-Hugoniot conditions are satisfied. Here, the flow first passes through the
outer saddle-type sonic point, then through a standing shock and finally through the
inner saddle-type sonic point. It was subsequently shown by MSC96 and
Ryu, Chakrabarti and Molteni (1997) that even when the parameters of the 
injected flow are outside this region oscillating shocks may form.
In the present paper we use Smoothed Particle Hydrodynamics (SPH) to do the 
simulations as in MSC96, but we use both halves of the flow 
to inject matter. The code has already been tested for its accuracy against 
theoretical solutions already (Chakrabarti \& Molteni, 1993, MSC96) and we do not repeat this here. 
We just want to recall that the `pseudo'-particles have been chosen 
to be `toroidal' in shape and they preserve angular momentum very 
accurately (Molteni, Ryu \& Chakrabarti, 1996). In MSC96, 
suggestion was made that {\it radial} oscillations of the shock may 
take place when the infall time in the post-shock region is comparable to the 
cooling time scale. Presently, we concentrate on the vast region of the parameter space 
where both the {\it radial and vertical} oscillations are present. 

In the context of galactic black holes one class of QPVs have been observed that are commonly 
known as Quasi-Periodic Oscillations or QPOs. Our solutions for QPVs are so generic that
we believe that QPVs should be observed even for super-massive black holes. 
Actually there are indications that QPOs in supermassive
black holes may have been observed (Halpern \& Marshall, 1996). Supermassive
black holes are believed to be fueled by low angular momentum matter generated from 
winds of nearby stars. So, a similar consideration 
as that for wind accreting galactic black holes will apply. As far as Comptonization 
is concerned, we use the following simple procedure: (a) We use an accretion rate ${\dot m}_b$
and run the code with bremsstrahlung cooling alone.  (b) We get the temperature 
of the electrons and the optical depth $\tau_b$ of the 'mean' size of the CENBOL region 
($X_s$ in Table 1). To get the corresponding  physical parameters (as if the Compton cooling were included
to begin with) (c) we assume that the injected soft
photons have energy $\sim 1$kev for the stellar mass black holes and $\sim 0.01$kev for  the 
massive black holes and using that (d) we compute the enhancement factor $\eta$ from DLC91 
for uniform spherical geometry for this injected photon and the electron temperature 
(assumed to be lower by a factor of $(m_p/m_e)^{1/2}$ than the proton temperature [e.g., Rees, 1984]) obtained
in each run. (e) We reduce the accretion 
rate ${\dot m}_b$ by this $\eta$ to get ${\dot m}_c={\dot m}_b/\eta$ and the optical
depth by $\tau_c=\tau_b/\eta$. For self-consistency, we ensure that $\tau_c$ thus obtained
is the same as that used while computing $\eta$ from DLC91. We also 
compute the Compton $y$ parameter to check if it has an acceptable value.
This procedure is re-formulates the problem of inclusion of the exact cooling
processes into the case where a simplified cooling effect
may be used. Since we are interested in the integrated energy loss due to cooling
and not the spectral properties, this procedure should give reasonable answer.

Table 1 lists the parameters for the Model runs 
we report in this paper. There are basically two groups of 
input. In Group A, we consider the black hole to be super-massive ($M=10^8M_\odot$) 
and in Group B, we consider a stellar mass black hole ($M=10M_\odot$). Different
runs are characterized by the injected matter density (gm/cc) at the outer boundary 
of the numerical grid given in the second column. In all the model runs, we choose 
the following parameters: the index $\gamma$ in the polytropic relation $P \propto \rho^{\gamma}$ 
($P$ is the isotropic pressure and $\rho$ is the gas density) is $5/3$, the outer 
boundary $r_{out}=50$, the specific angular momentum $\lambda=1.75$ (in units of $2GM/c$), 
injected radial velocity $\vartheta_r=0.126$, the sound velocity $a=0.04$ 
and the vertical height $h=15$. The vertical height at 
the outer boundary is so chosen that the flow remains in hydrostatic
equilibrium in the vertical direction at the point of injection 
$h \sim a r_{out}^{3/2}$. We measure all the 
distances in units of the Schwarzschild radius of the black hole 
$r_g=2GM_{BH}/c^2$, all the velocities in units of the velocity of 
light $c$ and all the masses in units of $M_{BH}$, the mass of the 
black hole. $c$ and $G$ are the velocity of light and the
gravitational constant respectively. Note that the 
marginally stable (lowest possible angular momentum with a stable Keplerian
orbit) angular momentum is $1.83$ in these units. 
In the third column we give the enhancement factor $\eta$ that we obtained using DLC91 for our run
and in the fourth column we provide ${\dot m}_c$ in units of ${\dot M_{Edd}}$,
the Eddington rate, the effective accretion rate
in the presence of Comptonization. In the fifth column we give the Compton $y$ parameters.
We assume soft injected photon of energy $h\nu_{inj} \sim 1$keV for galactic  black holes, and $0.01$keV 
for supermassive black holes in order to calculate $\eta$. In the sixth column,
we give the average observed location of the shock.
In the seventh column, we give the average number of `pseudo'-particles present in 
the disk in each run. In the eight column we present the QPV(QPO) frequency 
in Hz as obtained from PDS. In the ninth column, the inverse of the infall timescale
$t_i \sim X_s^{3/2}$ corrected for a constant compression ratio $R$ ($\sim 4$ for a
strong shock) is provided as an indication of the expected oscillation time scale.
In the tenth column, the ratio of the quantities in the eighth and ninth column
is given. In the eleventh column, the $rms$ amplitude in percentage is given,
we shall discuss in the next section. As in MSC96, 
we choose the Paczy\'nski-Wiita (1980) pseudo-Newtonian potential 
to describe the space-time around the Schwarzschild black hole. 
We believe that the result would remain unchanged if a true Schwarzschild
geometry were used, since the transonic properties do not change 
very much. However, if a Kerr geometry were used, we expect the locations of the
inner sonic points and shocks to move inward (Chakrabarti, 1996b) giving rise to higher QPO frequencies.


\centerline{ TABLE 1: Inputs and extracted parameters for the model runs}

\begin{center}
\begin{tabular}{lllllllllllll}
\hline
\hline
Run &$\rho_{inj}$ & $\eta$ & ${\dot m}_c$ & y & ${\bar X}_s$ & ${\bar N}$ & $\nu_{QPO}$ & $(Rt_i)^{-1}$ & $\frac{\nu_{QPO}}{(Rt_i)^{-1}}$ &  $R_{\it rms}$  \\
 & (gm/cc) &  &   &  & & & (Hz) & & & ($\%$) \\
\hline
\hline
A.1 & 0.95e-14 & 90 & 0.27 & 5.1& 16 & 22250 & 1.83e-6 & 3.9e-6 & 0.46  &  11.2 \\
&  & & &  & & & 9.15e-8 &  & &  \\
A.2 & 1.26e-14 & 142 & 0.22 & 6.6& 13 & 19250 &  2.06e-6 & 5.3e-6 & 0.39 & 11.2  \\
& &  & & &  &  & 5.34e-7 & & &   \\
A.3 & 3e-14 & 208&  0.37  & 8.3 & 6  & 12875 &  1.37e-6& 1.7e-5 & 0.23  &   11.2 \\
& &  & &  & & & 3.36e-7 & & &  \\
A.4 & 5e-14 & 331 & 0.39 & 9.4 & 5  & 10500 & 6.88e-7 & 2.2e-5 & 0.03 & 3.3 \\
&  & & & & & &  3.57e-7  & & &  \\
\hline
B.1  & 3.6e-10 & 1.8& 0.05 & 0.4& 15.8  & 21950  & 19.34  & 40 & 0.48 & 3.5\\
B.2   & 4.5e-10 & 1.5& 0.08 & 0.5 & 16.3 &22050  & 19.45 & 38.5  & 0.51 &  4.0 \\
B.3   & 4.5e-8 & 29.4 &0.39 & 11.8 & 24  & 33500  & 10.2& 21.3  & 0.48 &  14 \\
 & & &  & & & & 2.13 &  & & \\
B.4   & 4.5e-7 & 134 & 0.85 &43.4 &5 & 10825  &  8.32  & 223  & 0.04 & 3.7 \\
 & & & & &  & & 3.58 &  & &  \\
\hline
\hline
\end{tabular}
\end{center}

One of the important conditions of a standing shock
is that the thermal pressure ($P$) plus ram pressure ($\rho \vartheta_r^2$) should be 
continuous across the shock in the steady state, i.e.,
$$
P_- + \rho_- v_-^2 =P_+ + \rho_+ v_+^2 .
\eqno{(2)}
$$
Here, the subscripts $-$ and $+$ denote quantities before and after the shock 
and $v$ denotes the radial velocity of the flow as measured in local rotating frames.
When we increase the density of the gas, the cooling rate is increased. 
The thermal pressure is thus decreased, especially in the post-shock region,
causing the shock to move closer to the black hole.
In our simulations we generally see this 
trend in Group A runs (Table 1). Along with the shock location, there is a systematic 
decrease in frequency as well since the flow becomes cooler and the timescale of 
cooling itself is increased. 

In Group B Runs, within our parameter range, occasionally there were stronger 
shocks which oscillated radially, and the interaction between the outflow and the 
inflow produces vertical oscillations as well. Strong convection 
and turbulences were seen in the post-shock region which contributed to pushing the 
shock upstream. As a result, the shock location was not seen to change 
monotonically with accretion rate. The location is clearly determined by 
several effects: centrifugal force, thermal pressure, ram pressure, turbulent
pressure  and the cooling rate. The nature of the dependence has been discussed in earlier
work in detail (Chakrabarti, 1989; Chakrabarti \& Molteni, 1995; Molteni, Lanzafame \& 
Chakrabarti, 1994; MSC96) and we do not repeat this here. We will just recall that 
the shock location is independent of the location of the 
numerical boundary when the fundamental quantities such as specific angular 
momentum, accretion rate and specific energy at a given point (say, at the 
inner sonic point) are kept fixed.  Cases B.1 and B.2 have very small 
Compton parameters and negligible Comptonization. Cases B.3 and B.4 are 
presented for academic purposes, just to indicate the effects of excess cooling. 
The Compton parameters are very large in these two cases. 

For reference, we may add that the light crossing times
of the horizon for the two classes of black hole are $r_g/c \sim 10^3$s 
and $10^{-4}$s respectively. Since a steady shock (in the Schwarzschild geometry) may form anywhere between 
$X_s\sim 10$ to $100r_g$, the QPO time-period, if assumed to be 
related to the infall time $t \propto X_s^{3/2}$, may be close 
to $10^5$ to $10^8$s for super-massive black holes, and close to $10^{-2}$s 
to $10$s for stellar mass black holes. Since we chose the outer boundary of 
the simulation to be at $r_{out}=50$, the infall time is roughly $r_{out}^{3/2} 
\sim 350$ in units of $r_g/c$. In order to trust the results of our simulations, 
we ran the code for a duration several hundred times of this timescale (typically, $T_{run}\sim
50,000-60,000$ or more). For a black hole of $M=10M_\odot$, this time corresponds to only $5-6$s
of real time. 

\section{Results}

First, we show examples of shock formation which exhibit radial and/or 
vertical oscillations.  Fig. 1(a-c) shows the locations of the SPH 
particles (dots) along with the velocity vectors (arrows) of 
every fifth particle for clarity. This case corresponds to 
case A.2. The matter distribution is shown at three different times (in
units of $r_g/c$) illustrating the vertical and radial motions of 
the shock which oscillates  between $11$ to $15$ Schwarzschild radii, mostly 
staying close to $r_s\sim 13r_g$. Some matter can be seen bouncing back
from the centrifugal barrier (boundary between the empty region along the vertical axis and the disk/jet 
matter) near the axis and forming giant vortices which interact with the 
inflow. These vortices push the post-shock region ($r\lsim 15r_g$) 
alternately into the upper and the lower halves.

\begin {figure}
\vbox{
\vskip -3.0cm
\hskip 0.0cm
\centerline{
\psfig{figure=f1.eps,height=13truecm,width=13truecm}}}
\vspace{0.0cm}
\caption[] { Snapshots of simulations of accretion disks around a $10^8M_\odot$ black hole
by Smoothed Particle Hydrodynamics. The dots are particle locations and arrows are
drawn for every fifth particle for clarity. Time (in units of $r_g/c$) is marked
in each box. Note the vertical as well as radial oscillation of the accretion shock wave
located at $\sim 13r_g$.}
\end{figure}

In Fig. 2(a-d) we show the variation of the light curve in runs A.1-A.4 (marked 
by a-d in the Figure). The light curves are  the variations of the
luminosity of bremsstrahlung radiation (in arbitrary units) for all the
matter within $r=50$ flow as a function of time (in seconds). 
As the density is increased, the average thermal energy gradually goes down
due to the presence of enhanced cooling and thus the luminosity also
goes down. The average location of the shock decreases (MSC96, Table 1).
The average number of pseudo-particles (Table 1) as well as the variation of the 
number of these particles also go down. As a result, with  
the increase in cooling loss, the light curves are 
found to be less noisy, and the amplitude of fluctuation is found to be lower. 

\begin {figure}
\vbox{
\vskip -0.0cm
\hskip 0.0cm
\centerline{
\psfig{figure=f2.eps,height=13truecm,width=13truecm}}}
\vspace{-2.0cm}
\noindent {\small {\bf Fig. 2(a-d)}}:
Total bremsstrahlung luminosity of the accretion flow
(in arbitrary units) as a function of time (in seconds).
The injected densities are (a) $\rho_{inj}=0.95\times 10^{-14}$gm/s,
(b) $1.265 \times 10^{-14}$gm/s, (c) $3 \times 10^{-14}$gm/s and
(d) $5 \times 10^{-14}$gm/s respectively. The disk becomes
cooler with increasing accretion rate due to bremsstrahlung
energy loss.
\end{figure}

In Fig. 3(a-d) we show the results of the Fast Fourier Transform of these light curves using FTOOLS
provided by NASA, the same software package that is used to 
analyze observational results. We find the remarkable result of familiar
power density spectra (PDS) with a QPV near the break frequency.
The QPV frequencies, indicated by the arrows, are presented in Table 1. 
We find that, perhaps due to the
presence of both the vertical and radial oscillations, there are 
two QPVs; the frequency of the stronger one gradually increases with 
the cooling, while that of the weaker one gradually decreases. 
The occurrence of QPV frequencies  at or near the break frequency of the PDS
can be understood by the following: according
to our solution, QPVs occur due to oscillations of the shock waves that separate
two `phases' of the flow --- the pre-shock flow is supersonic and weakly radiating,
while the post-shock flow is subsonic and strongly radiating. The former 
region of the flow produces the `flat-top' region in the PDS, while the latter region produces
a PDS with a slope $\sim -2$. The rms amplitude of the flat top
region up to the break frequency is shown in Table 1. It has generally similar value at around $11.2$\% except 
in (d) where it is only around $3.3$\%. These rms values are similar to the observed results from black hole 
candidates as well. In Fig. 3(a-b), where the cooling is weaker, the QPVs are less prominent
and we had to search the locally significant peaks of the PDS for them. As the cooling becomes stronger, 
the peak becomes stronger. 

\begin {figure}
\vbox{
\vskip -0.0cm
\hskip 0.0cm
\centerline{
\psfig{figure=f3.eps,height=13truecm,width=13truecm}}}
\vspace{-2.0cm}
\noindent {\small {\bf Fig. 3(a-d)}}: Power Density Spectra (PDS) of the four cases shown in Fig. 2.
Quasi Periodic Variation frequencies (denoted by $\nu_{QPO}$) can be
seen with timescales of hours to weeks. QPV peaks
(marked with arrows) are located near break frequencies with flat top
behavior at low frequency and power-law behavior
at high frequency. See Table 1 for other properties.
\end{figure}

In Fig. 4(a-d), we present similar results as in Fig. 2 for stellar mass black holes
(runs B.1 - B.4 of Table 1). Here, too, the disk becomes cooler with increasing density 
and the light curve changes its character according to whether the shocks are strong or not. 
In the beginning of the simulation some transient behaviour is seen till the 
disk settles down within a second in real time.
In Fig. 4(c), the parameters are such that a strong shock forms and it heats up the
disk while for the parameters in Fig. 4(d) a very weak shock forms. Shock oscillation
produces a large amplitude noise in the B.3 case (Fig. 4c). In Fig. 5(a-d),
the PDS of these four runs are shown. The transient region has been excluded 
while doing the FFT. As in Figs. 3(a-d), here too QPVs are
observed. For small black holes QPVs occur at a few Hz and the QPV is located at or 
near the break frequency. Thus the characteristics are similar to those of QPOs observed from 
black hole candidates. Radial and vertical oscillations often produce 
multiple peaks in this case as well. Other details can be seen in Table 1.
In Fig. 5c the peak is not very strong, perhaps due to the large amplitude noise
that is produced in this case (See, Fig. 4c). 

\begin {figure}
\vbox{
\vskip -0.0cm
\hskip 0.0cm
\centerline{
\psfig{figure=f4.eps,height=13truecm,width=13truecm}}}
\vspace{-2.0cm}
\noindent {\small {\bf Fig. 4(a-d)}}: Total bremsstrahlung luminosity of the accretion flow
(in arbitrary units) as a function of time (in seconds).
The injected densities are (a) $\rho_{inj}=3.6 \times 10^{-10}$gm/s,
(b) $4.5\times 10^{-10}$gm/s, (c) $4.5 \times 10^{-8}$gm/s and
(d) $4.5 \times 10^{-7}$gm/s respectively. The disk becomes
cooler with increasing accretion rate due to bremsstrahlung
energy loss.
\end{figure}

Occasionally, peaks at a very high frequency could be seen, but these peaks are
found to be transient. If the entire light curve is broken into smaller
pieces, weak peaks at frequencies $\sim 100-300$Hz can be seen, though they
do not persist. One example of this is shown in Fig. 6(a) where the PDS of the 
average thermal energy up to $T=36000$ is plotted for the case B.4. In Fig. 6(b) and (c) 
we show the flow pattern at two times which are separated by only $18$ units
($0.0018$s). The shock located at $\sim 6r_g$ shows two distinctly different
shapes at these two times, bent up and bent down. Oscillation of a density enhancement at $\sim 6r_g$ 
would have a time period of $\sim 6^{3/2} \times 10^{-4} R \sim 0.0015 R $s, where $R$, the compression
ratio, is $\sim 2 $ for a weak shock. The corresponding frequency is $333$Hz. The 
weak peak at $306\pm 7.93$Hz (Fig. 6a) may be due to this oscillation.
When data for much longer time is analyzed, this weak peak disappears completely.

\begin {figure}
\vbox{
\vskip -0.0cm
\hskip 0.0cm
\centerline{
\psfig{figure=f5.eps,height=13truecm,width=13truecm}}}
\vspace{-2.0cm}
\noindent {\small {\bf Fig. 5(a-d)}}: 
Power Density Spectra (PDS) of the four cases shown in Fig. 4.
Quasi Periodic Variation frequencies (denoted by $\nu_{QPO}$) can be
seen in the range of $2-10$Hz and also
close to several hundred Hz. QPV peaks (marked with arrows)
are located near break frequencies with flat top
behavior at low frequency and power-law behavior
at high frequency. See Table 1 for other properties.
\end{figure}


\begin {figure}
\vbox{
\vskip -0.0cm
\hskip 0.0cm
\centerline{
\psfig{figure=f6.eps,height=13truecm,width=13truecm}}}
\vspace{-0.0cm}
\noindent {\small {\bf Fig. 6(a-c)}}: 
(a) Power density spectrum of the average thermal energy for case B.4 with time up to
3.6s. Snapshots of simulations at closely separated times (marked) are
shown in (b-c). Only the inner $10$ Schwarzschild radii are shown.
The dots are particle locations and arrows are
drawn for every fifth particle for clarity. Time (in units of $r_g/c$) is marked
in each box. The oscillation of this region causes QPV at $\sim 300$Hz [marked with an arrow
in (a)].
\end{figure}

\section{The possible cause of shock behaviour and the resulting QPV}

A shock remains steady, or `standing' if Eq. (2) along with the 
specific energy and angular momentum flux remain continuous across it.
These conditions are known as the Rankine-Hugoniot 
conditions (Landau \& Lifshitz, 1959). However, in the presence 
of an efficient energy extraction mechanism from the disk, whether through thermal/non-thermal radiation 
or through higher entropy flows along the shocked winds and outflows, the shock  should
move closer to the black hole as the pressure in the post-shock region ($P_+$) drops.
This is what is seen  (see, Table 1), except in case B.3 where the turbulence was
too strong to push the shock outward. 

What about the oscillation time scale?  We note that the cooling time 
scale $ t_c\sim {\cal E}/{\dot {\cal E}}$ (where ${\cal E}$ is the 
specific energy and a `dot' indicates its time derivative) competes 
with the infall time scale $t_i \sim X_s/v$. MSC96 argued that 
when these two time scales agree roughly, i.e., roughly when the resonance takes place,
the shocks may start oscillating. The infall time $t_i$ is not easy to compute,
however. CM00 argued that the infall velocity need not be free-fall 
velocity $v\propto 1/x^{1/2}$ as the post-shock flow is reduced at the shock by at least 
a factor of $R$, the compression ratio. Furthermore, the velocity is 
highly modified due to turbulent pressure and angular motion and could be very slowly varying
in between the shock and the inner sonic point (CM00; Chakrabarti, Ryu \& Samanta, in preparation). 
A third and perhaps the most important parameter in predicting the time scale
is the rate of outflow from the post-shock region and the extent to which 
the total mass of the disk oscillates. These two cooling mechanisms, namely, through radiation loss
and through outflow, are acting in opposite directions -- enhanced thermal cooling
reduces the pressure gradient force, thereby reducing the wind rate. 
The QPV timescale is expected to be close to  
$t_{QPV}=\nu_{QPV}^{-1}\sim t_c \sim Rt_i$. In Table 1, $\nu_{QPV}/(Rt_i)^{-1}$
is presented. This ratio seems to be closer to $0.5$ (justifying a {\it near} resonance condition)
except when the cooling is very strong (runs A.4 and B.4) when the
ratio drastically goes down. In these latter types of runs the shock moves very close to the
black hole and becomes very weak. The QPV could be due to the domination of the
interaction of the winds with the inflow. The two QPV frequencies
often seen could be due to the two types of cooling.

\section{Concluding remarks}

In this paper, we presented the discovery of Quasi-Periodic 
Variations of emitted radiation  through extensive time-dependent 
numerical simulations. In earlier studies, such as in MSC96 and CM00, it was
suggested that shock oscillations might contribute to the observed QPOs. 
In the present work, we studied shock oscillations in low-angular-momentum 
axisymmetric flows around stellar and super-massive black holes for a range of accretion rates.
Unlike MSC96, we used fully two dimensional inflow (upper and lower halves)
and demonstrated that the shocks participate in vertical and radial oscillations. 
Using PDS of light curves obtained from time dependent numerical 
simulations of sub-Keplerian flows, we demonstrated that the  QPVs 
could be the cause of the observed QPOs.
We showed that the computed light curve produces power density spectra with 
characteristic QPVs located at or near the break frequency which separates the flat-top 
type and power-law type PDSs. For stellar black holes, we see QPVs at 
reasonable frequencies ($2-20$Hz). Very often high frequency QPVs 
(at around $100-300$Hz) could be observed. We demonstrated that
this could be due to oscillation at the very inner edge where the flow becomes supersonic.
For Active Galaxies and Quasars, QPVs at very low frequencies ($\sim 10^{-7}-10^{-6}$ Hz)
are observed. Their properties are similar to those of QPOs observed in these candidates.

Based on our simulations we propose that the QPOs observed in black hole 
candidates are genuinely related to the accretion shock oscillations. 
Variation of QPO frequencies from time to time in the same object are thus
manifestations of the variation of the Keplerian and sub-Keplerian accretion rates. 
In our work, we assumed a simple  Paczy\'nski-Wiita potential 
to study the shock properties around Schwarzschild black holes. Based on our experience
with transonic flows, we are confident that none of the properties would be affected 
if the simulations were carried out in true Schwarzschild geometry. However, if a rotating
black hole is used, the shocks will form much closer to the black hole (Chakrabarti, 1996b), resulting in
higher QPV frequencies. 

SKC thanks the support of Indian Space Research Organization for a RESPOND project. 
Discussions with Prof. A.R. Rao (TIFR) and Mr. A. Nandi (SNBNCBS) are also acknowledged. 

\begin{thebibliography}{}
\def\ref#1\par{\parshape=2 0in 14.5cm 1cm 13.5cm {#1} \par}
\parskip=0pt
\parindent=0pt

\bibitem[]{}
Chakrabarti S.K. 1989, ApJ, 347, 365

\bibitem[]{}
Chakrabarti S.K. 1996a, ApJ, 464, 664

\bibitem[]{}
Chakrabarti S.K. 1996b, MNRAS, 283, 325

\bibitem[]{}
Chakrabarti S.K. 1999, A\&A 351, 185

\bibitem[]{}
Chakrabarti S.K., Manickam, S.G. 2000, ApJ 531, L41 (CM00)

\bibitem[]{}
Chakrabarti S.K., Molteni, D.G. 1995, MNRAS, 272, 80

\bibitem[]{} 
Chakrabarti, S.K.,  Nandi, A., Manickam, S. G.,  Mandal, S.  Rao, A.R. 2002, ApJ 579, L21

\bibitem[]{}
Chanmugam, G., Langer, S.H., Shaviv, G. 1985, ApJ, 288, 428

\bibitem[]{}
Dermer, C.D., Liang, E.P., Canfield, E. 1991, ApJ, 369, 410

\bibitem[]{}
Dhawan, V., Mirabel, I. F., Rodriguez, L. F. 2000, ApJ 543 373

\bibitem[]{}
Halpern, J.P., Marshall, H.L. 1996,  ApJ, 464, 760

\bibitem[]{}
Landau, L.D., Lifshitz, E.D. 1959, Fluid Mechanics (New York:Pergamon)

\bibitem[]{}
Langer, S.H., Chanmugam, G., Shaviv, G. 1981, ApJ, 245, L23

\bibitem[]{} 
Molteni, D., Acharyya, K., Kuznetsov, O., Bisikalo, D.,  Chakrabarti, S.\ K. 2001, ApJL, 563, L57

\bibitem[]{} 
Molteni, D., Sponholz, H.,  Chakrabarti, S.\ K. 1996, ApJ, 457, 805  (MSC96)

\bibitem[]{} 
Molteni, D., Ryu, D., Chakrabarti, S.\ K. 1996, ApJ, 470, 460

\bibitem[]{}
Paczy\'nski, B., Wiita, P.J., 1980, A\&A, 88, 23 

\bibitem[]{}
Rao, A.R., Naik, S., Vadawale, S.V., Chakrabarti, S.K. 2000, A\&A 360, L25

\bibitem[]{}
Rees, M.J. 1984, Ann. Rev. Astron. Astrophys., 22, 471

\bibitem[]{}
Rybicki, G., Lightman, A.P., 1979, Radiative Processes in Astrophysics (New York:Wiley Interscience)

\bibitem[]{}
Ryu, D., Chakrabarti, S.K., Molteni, D. 1997, ApJ, 474, 378 (RCM97)

\bibitem[]{}
Shaham, J. 1987 in The Origin and Evolution of Neutron Stars, D.J. Helfand and J.-H. Huang (eds.), 347

\bibitem[]{}
Smith, D. M., Heindl, W. A., Swank, J. H., 2002, ApJ, 569, 362

\end{thebibliography}{}

\end{document}


psfig.tex
% Psfig/TeX 
\def\PsfigVersion{1.9}
% dvips version
%
% All psfig/tex software, documentation, and related files
% in this distribution of psfig/tex are 
% Copyright 1987, 1988, 1991 Trevor J. Darrell
%
% Permission is granted for use and non-profit distribution of psfig/tex 
% providing that this notice is clearly maintained. The right to
% distribute any portion of psfig/tex for profit or as part of any commercial
% product is specifically reserved for the author(s) of that portion.
%
% *** Feel free to make local modifications of psfig as you wish,
% *** but DO NOT post any changed or modified versions of ``psfig''
% *** directly to the net. Send them to me and I'll try to incorporate
% *** them into future versions. If you want to take the psfig code 
% *** and make a new program (subject to the copyright above), distribute it, 
% *** (and maintain it) that's fine, just don't call it psfig.
%
% Bugs and improvements to trevor@media.mit.edu.
%
% Thanks to Greg Hager (GDH) and Ned Batchelder for their contributions
% to the original version of this project.
%
% Modified by J. Daniel Smith on 9 October 1990 to accept the
% %%BoundingBox: comment with or without a space after the colon.  Stole
% file reading code from Tom Rokicki's EPSF.TEX file (see below).
%
% More modifications by J. Daniel Smith on 29 March 1991 to allow the
% the included PostScript figure to be rotated.  The amount of
% rotation is specified by the "angle=" parameter of the \psfig command.
%
% Modified by Robert Russell on June 25, 1991 to allow users to specify
% .ps filenames which don't yet exist, provided they explicitly provide
% boundingbox information via the \psfig command. Note: This will only work
% if the "file=" parameter follows all four "bb???=" parameters in the
% command. This is due to the order in which psfig interprets these params.
%
%  3 Jul 1991	JDS	check if file already read in once
%  4 Sep 1991	JDS	fixed incorrect computation of rotated
%			bounding box
% 25 Sep 1991	GVR	expanded synopsis of \psfig
% 14 Oct 1991	JDS	\fbox code from LaTeX so \psdraft works with TeX
%			changed \typeout to \ps@typeout
% 17 Oct 1991	JDS	added \psscalefirst and \psrotatefirst
%

% From: gvr@cs.brown.edu (George V. Reilly)
%
% \psdraft	draws an outline box, but doesn't include the figure
%		in the DVI file.  Useful for previewing.
%
% \psfull	includes the figure in the DVI file (default).
%
% \psscalefirst width= or height= specifies the size of the figure
% 		before rotation.
% \psrotatefirst (default) width= or height= specifies the size of the
% 		 figure after rotation.  Asymetric figures will
% 		 appear to shrink.
%
% \psfigurepath#1	sets the path to search for the figure
%
% \psfig
% usage: \psfig{file=, figure=, height=, width=,
%			bbllx=, bblly=, bburx=, bbury=,
%			rheight=, rwidth=, clip=, angle=, silent=}
%
%	"file" is the filename.  If no path name is specified and the
%		file is not found in the current directory,
%		it will be looked for in directory \psfigurepath.
%	"figure" is a synonym for "file".
%	By default, the width and height of the figure are taken from
%		the BoundingBox of the figure.
%	If "width" is specified, the figure is scaled so that it has
%		the specified width.  Its height changes proportionately.
%	If "height" is specified, the figure is scaled so that it has
%		the specified height.  Its width changes proportionately.
%	If both "width" and "height" are specified, the figure is scaled
%		anamorphically.
%	"bbllx", "bblly", "bburx", and "bbury" control the PostScript
%		BoundingBox.  If these four values are specified
%               *before* the "file" option, the PSFIG will not try to
%               open the PostScript file.
%	"rheight" and "rwidth" are the reserved height and width
%		of the figure, i.e., how big TeX actually thinks
%		the figure is.  They default to "width" and "height".
%	The "clip" option ensures that no portion of the figure will
%		appear outside its BoundingBox.  "clip=" is a switch and
%		takes no value, but the `=' must be present.
%	The "angle" option specifies the angle of rotation (degrees, ccw).
%	The "silent" option makes \psfig work silently.
%

% check to see if macros already loaded in (maybe some other file says
% "\input psfig") ...
\ifx\undefined\psfig\else\endinput\fi

%
% from a suggestion by eijkhout@csrd.uiuc.edu to allow
% loading as a style file. Changed to avoid problems
% with amstex per suggestion by jbence@math.ucla.edu

\let\LaTeXAtSign=\@
\let\@=\relax
\edef\psfigRestoreAt{\catcode`\@=\number\catcode`@\relax}
%\edef\psfigRestoreAt{\catcode`@=\number\catcode`@\relax}
\catcode`\@=11\relax
\newwrite\@unused
\def\ps@typeout#1{{\let\protect\string\immediate\write\@unused{#1}}}
\ps@typeout{psfig/tex \PsfigVersion}

%% Here's how you define your figure path.  Should be set up with null
%% default and a user useable definition.

\def\figurepath{./}
\def\psfigurepath#1{\edef\figurepath{#1}}

%
% @psdo control structure -- similar to Latex @for.
% I redefined these with different names so that psfig can
% be used with TeX as well as LaTeX, and so that it will not 
% be vunerable to future changes in LaTeX's internal
% control structure,
%
\def\@nnil{\@nil}
\def\@empty{}
\def\@psdonoop#1\@@#2#3{}
\def\@psdo#1:=#2\do#3{\edef\@psdotmp{#2}\ifx\@psdotmp\@empty \else
    \expandafter\@psdoloop#2,\@nil,\@nil\@@#1{#3}\fi}
\def\@psdoloop#1,#2,#3\@@#4#5{\def#4{#1}\ifx #4\@nnil \else
       #5\def#4{#2}\ifx #4\@nnil \else#5\@ipsdoloop #3\@@#4{#5}\fi\fi}
\def\@ipsdoloop#1,#2\@@#3#4{\def#3{#1}\ifx #3\@nnil 
       \let\@nextwhile=\@psdonoop \else
      #4\relax\let\@nextwhile=\@ipsdoloop\fi\@nextwhile#2\@@#3{#4}}
\def\@tpsdo#1:=#2\do#3{\xdef\@psdotmp{#2}\ifx\@psdotmp\@empty \else
    \@tpsdoloop#2\@nil\@nil\@@#1{#3}\fi}
\def\@tpsdoloop#1#2\@@#3#4{\def#3{#1}\ifx #3\@nnil 
       \let\@nextwhile=\@psdonoop \else
      #4\relax\let\@nextwhile=\@tpsdoloop\fi\@nextwhile#2\@@#3{#4}}
% 
% \fbox is defined in latex.tex; so if \fbox is undefined, assume that
% we are not in LaTeX.
% Perhaps this could be done better???
\ifx\undefined\fbox
% \fbox code from modified slightly from LaTeX
\newdimen\fboxrule
\newdimen\fboxsep
\newdimen\ps@tempdima
\newbox\ps@tempboxa
\fboxsep = 3pt
\fboxrule = .4pt
\long\def\fbox#1{\leavevmode\setbox\ps@tempboxa\hbox{#1}\ps@tempdima\fboxrule
    \advance\ps@tempdima \fboxsep \advance\ps@tempdima \dp\ps@tempboxa
   \hbox{\lower \ps@tempdima\hbox
  {\vbox{\hrule height \fboxrule
          \hbox{\vrule width \fboxrule \hskip\fboxsep
          \vbox{\vskip\fboxsep \box\ps@tempboxa\vskip\fboxsep}\hskip 
                 \fboxsep\vrule width \fboxrule}
                 \hrule height \fboxrule}}}}
\fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% file reading stuff from epsf.tex
%   EPSF.TEX macro file:
%   Written by Tomas Rokicki of Radical Eye Software, 29 Mar 1989.
%   Revised by Don Knuth, 3 Jan 1990.
%   Revised by Tomas Rokicki to accept bounding boxes with no
%      space after the colon, 18 Jul 1990.
%   Portions modified/removed for use in PSFIG package by
%      J. Daniel Smith, 9 October 1990.
%
\newread\ps@stream
\newif\ifnot@eof       % continue looking for the bounding box?
\newif\if@noisy        % report what you're making?
\newif\if@atend        % %%BoundingBox: has (at end) specification
\newif\if@psfile       % does this look like a PostScript file?
%
% PostScript files should start with `%!'
%
{\catcode`\%=12\global\gdef\epsf@start{%!}}
\def\epsf@PS{PS}
%
\def\epsf@getbb#1{%
%
%   The first thing we need to do is to open the
%   PostScript file, if possible.
%
\openin\ps@stream=#1
\ifeof\ps@stream\ps@typeout{Error, File #1 not found}\else
%
%   Okay, we got it. Now we'll scan lines until we find one that doesn't
%   start with %. We're looking for the bounding box comment.
%
   {\not@eoftrue \chardef\other=12
    \def\do##1{\catcode`##1=\other}\dospecials \catcode`\ =10
    \loop
       \if@psfile
	  \read\ps@stream to \epsf@fileline
       \else{
	  \obeyspaces
          \read\ps@stream to \epsf@tmp\global\let\epsf@fileline\epsf@tmp}
       \fi
       \ifeof\ps@stream\not@eoffalse\else
%
%   Check the first line for `%!'.  Issue a warning message if its not
%   there, since the file might not be a PostScript file.
%
       \if@psfile\else
       \expandafter\epsf@test\epsf@fileline:. \\%
       \fi
%
%   We check to see if the first character is a % sign;
%   if so, we look further and stop only if the line begins with
%   `%%BoundingBox:' and the `(atend)' specification was not found.
%   That is, the only way to stop is when the end of file is reached,
%   or a `%%BoundingBox: llx lly urx ury' line is found.
%
          \expandafter\epsf@aux\epsf@fileline:. \\%
       \fi
   \ifnot@eof\repeat
   }\closein\ps@stream\fi}%
%
% This tests if the file we are reading looks like a PostScript file.
%
\long\def\epsf@test#1#2#3:#4\\{\def\epsf@testit{#1#2}
			\ifx\epsf@testit\epsf@start\else
\ps@typeout{Warning! File does not start with `\epsf@start'.  It may not be a PostScript file.}
			\fi
			\@psfiletrue} % don't test after 1st line
%
%   We still need to define the tricky \epsf@aux macro. This requires
%   a couple of magic constants for comparison purposes.
%
{\catcode`\%=12\global\let\epsf@percent=%\global\def\epsf@bblit{%BoundingBox}}
%
%
%   So we're ready to check for `%BoundingBox:' and to grab the
%   values if they are found.  We continue searching if `(at end)'
%   was found after the `%BoundingBox:'.
%
\long\def\epsf@aux#1#2:#3\\{\ifx#1\epsf@percent
   \def\epsf@testit{#2}\ifx\epsf@testit\epsf@bblit
	\@atendfalse
        \epsf@atend #3 . \\%
	\if@atend	
	   \if@verbose{
		\ps@typeout{psfig: found `(atend)'; continuing search}
	   }\fi
        \else
        \epsf@grab #3 . . . \\%
        \not@eoffalse
        \global\no@bbfalse
        \fi
   \fi\fi}%
%
%   Here we grab the values and stuff them in the appropriate definitions.
%
\def\epsf@grab #1 #2 #3 #4 #5\\{%
   \global\def\epsf@llx{#1}\ifx\epsf@llx\empty
      \epsf@grab #2 #3 #4 #5 .\\\else
   \global\def\epsf@lly{#2}%
   \global\def\epsf@urx{#3}\global\def\epsf@ury{#4}\fi}%
%
% Determine if the stuff following the %%BoundingBox is `(atend)'
% J. Daniel Smith.  Copied from \epsf@grab above.
%
\def\epsf@atendlit{(atend)} 
\def\epsf@atend #1 #2 #3\\{%
   \def\epsf@tmp{#1}\ifx\epsf@tmp\empty
      \epsf@atend #2 #3 .\\\else
   \ifx\epsf@tmp\epsf@atendlit\@atendtrue\fi\fi}


% End of file reading stuff from epsf.tex
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% trigonometry stuff from "trig.tex"
\chardef\psletter = 11 % won't conflict with \begin{letter} now...
\chardef\other = 12

\newif \ifdebug %%% turn me on to see TeX hard at work ...
\newif\ifc@mpute %%% don't need to compute some values
\c@mputetrue % but assume that we do

\let\then = \relax
\def\r@dian{pt }
\let\r@dians = \r@dian
\let\dimensionless@nit = \r@dian
\let\dimensionless@nits = \dimensionless@nit
\def\internal@nit{sp }
\let\internal@nits = \internal@nit
\newif\ifstillc@nverging
\def \Mess@ge #1{\ifdebug \then \message {#1} \fi}

{ %%% Things that need abnormal catcodes %%%
	\catcode `\@ = \psletter
	\gdef \nodimen {\expandafter \n@dimen \the \dimen}
	\gdef \term #1 #2 #3%
	       {\edef \t@ {\the #1}%%% freeze parameter 1 (count, by value)
		\edef \t@@ {\expandafter \n@dimen \the #2\r@dian}%
				   %%% freeze parameter 2 (dimen, by value)
		\t@rm {\t@} {\t@@} {#3}%
	       }
	\gdef \t@rm #1 #2 #3%
	       {{%
		\count 0 = 0
		\dimen 0 = 1 \dimensionless@nit
		\dimen 2 = #2\relax
		\Mess@ge {Calculating term #1 of \nodimen 2}%
		\loop
		\ifnum	\count 0 < #1
		\then	\advance \count 0 by 1
			\Mess@ge {Iteration \the \count 0 \space}%
			\Multiply \dimen 0 by {\dimen 2}%
			\Mess@ge {After multiplication, term = \nodimen 0}%
			\Divide \dimen 0 by {\count 0}%
			\Mess@ge {After division, term = \nodimen 0}%
		\repeat
		\Mess@ge {Final value for term #1 of 
				\nodimen 2 \space is \nodimen 0}%
		\xdef \Term {#3 = \nodimen 0 \r@dians}%
		\aftergroup \Term
	       }}
	\catcode `\p = \other
	\catcode `\t = \other
	\gdef \n@dimen #1pt{#1} %%% throw away the ``pt''
}

\def \Divide #1by #2{\divide #1 by #2} %%% just a synonym

\def \Multiply #1by #2%%% allows division of a dimen by a dimen
       {{%%% should really freeze parameter 2 (dimen, passed by value)
	\count 0 = #1\relax
	\count 2 = #2\relax
	\count 4 = 65536
	\Mess@ge {Before scaling, count 0 = \the \count 0 \space and
			count 2 = \the \count 2}%
	\ifnum	\count 0 > 32767 %%% do our best to avoid overflow
	\then	\divide \count 0 by 4
		\divide \count 4 by 4
	\else	\ifnum	\count 0 < -32767
		\then	\divide \count 0 by 4
			\divide \count 4 by 4
		\else
		\fi
	\fi
	\ifnum	\count 2 > 32767 %%% while retaining reasonable accuracy
	\then	\divide \count 2 by 4
		\divide \count 4 by 4
	\else	\ifnum	\count 2 < -32767
		\then	\divide \count 2 by 4
			\divide \count 4 by 4
		\else
		\fi
	\fi
	\multiply \count 0 by \count 2
	\divide \count 0 by \count 4
	\xdef \product {#1 = \the \count 0 \internal@nits}%
	\aftergroup \product
       }}

\def\r@duce{\ifdim\dimen0 > 90\r@dian \then   % sin(x+90) = sin(180-x)
		\multiply\dimen0 by -1
		\advance\dimen0 by 180\r@dian
		\r@duce
	    \else \ifdim\dimen0 < -90\r@dian \then  % sin(-x) = sin(360+x)
		\advance\dimen0 by 360\r@dian
		\r@duce
		\fi
	    \fi}

\def\Sine#1%
       {{%
	\dimen 0 = #1 \r@dian
	\r@duce
	\ifdim\dimen0 = -90\r@dian \then
	   \dimen4 = -1\r@dian
	   \c@mputefalse
	\fi
	\ifdim\dimen0 = 90\r@dian \then
	   \dimen4 = 1\r@dian
	   \c@mputefalse
	\fi
	\ifdim\dimen0 = 0\r@dian \then
	   \dimen4 = 0\r@dian
	   \c@mputefalse
	\fi
%
	\ifc@mpute \then
        	% convert degrees to radians
		\divide\dimen0 by 180
		\dimen0=3.141592654\dimen0
%
		\dimen 2 = 3.1415926535897963\r@dian %%% a well-known constant
		\divide\dimen 2 by 2 %%% we only deal with -pi/2 : pi/2
		\Mess@ge {Sin: calculating Sin of \nodimen 0}%
		\count 0 = 1 %%% see power-series expansion for sine
		\dimen 2 = 1 \r@dian %%% ditto
		\dimen 4 = 0 \r@dian %%% ditto
		\loop
			\ifnum	\dimen 2 = 0 %%% then we've done
			\then	\stillc@nvergingfalse 
			\else	\stillc@nvergingtrue
			\fi
			\ifstillc@nverging %%% then calculate next term
			\then	\term {\count 0} {\dimen 0} {\dimen 2}%
				\advance \count 0 by 2
				\count 2 = \count 0
				\divide \count 2 by 2
				\ifodd	\count 2 %%% signs alternate
				\then	\advance \dimen 4 by \dimen 2
				\else	\advance \dimen 4 by -\dimen 2
				\fi
		\repeat
	\fi		
			\xdef \sine {\nodimen 4}%
       }}

% Now the Cosine can be calculated easily by calling \Sine
\def\Cosine#1{\ifx\sine\UnDefined\edef\Savesine{\relax}\else
		             \edef\Savesine{\sine}\fi
	{\dimen0=#1\r@dian\advance\dimen0 by 90\r@dian
	 \Sine{\nodimen 0}
	 \xdef\cosine{\sine}
	 \xdef\sine{\Savesine}}}	      
% end of trig stuff
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\def\psdraft{
	\def\@psdraft{0}
	%\ps@typeout{draft level now is \@psdraft \space . }
}
\def\psfull{
	\def\@psdraft{100}
	%\ps@typeout{draft level now is \@psdraft \space . }
}

\psfull

\newif\if@scalefirst
\def\psscalefirst{\@scalefirsttrue}
\def\psrotatefirst{\@scalefirstfalse}
\psrotatefirst

\newif\if@draftbox
\def\psnodraftbox{
	\@draftboxfalse
}
\def\psdraftbox{
	\@draftboxtrue
}
\@draftboxtrue

\newif\if@prologfile
\newif\if@postlogfile
\def\pssilent{
	\@noisyfalse
}
\def\psnoisy{
	\@noisytrue
}
\psnoisy
%%% These are for the option list.
%%% A specification of the form a = b maps to calling \@p@@sa{b}
\newif\if@bbllx
\newif\if@bblly
\newif\if@bburx
\newif\if@bbury
\newif\if@height
\newif\if@width
\newif\if@rheight
\newif\if@rwidth
\newif\if@angle
\newif\if@clip
\newif\if@verbose
\def\@p@@sclip#1{\@cliptrue}


\newif\if@decmpr

%%% GDH 7/26/87 -- changed so that it first looks in the local directory,
%%% then in a specified global directory for the ps file.
%%% RPR 6/25/91 -- changed so that it defaults to user-supplied name if
%%% boundingbox info is specified, assuming graphic will be created by
%%% print time.
%%% TJD 10/19/91 -- added bbfile vs. file distinction, and @decmpr flag

\def\@p@@sfigure#1{\def\@p@sfile{null}\def\@p@sbbfile{null}
	        \openin1=#1.bb
		\ifeof1\closein1
	        	\openin1=\figurepath#1.bb
			\ifeof1\closein1
			        \openin1=#1
				\ifeof1\closein1%
				       \openin1=\figurepath#1
					\ifeof1
					   \ps@typeout{Error, File #1 not found}
						\if@bbllx\if@bblly
				   		\if@bburx\if@bbury
			      				\def\@p@sfile{#1}%
			      				\def\@p@sbbfile{#1}%
							\@decmprfalse
				  	   	\fi\fi\fi\fi
					\else\closein1
				    		\def\@p@sfile{\figurepath#1}%
				    		\def\@p@sbbfile{\figurepath#1}%
						\@decmprfalse
	                       		\fi%
			 	\else\closein1%
					\def\@p@sfile{#1}
					\def\@p@sbbfile{#1}
					\@decmprfalse
			 	\fi
			\else
				\def\@p@sfile{\figurepath#1}
				\def\@p@sbbfile{\figurepath#1.bb}
				\@decmprtrue
			\fi
		\else
			\def\@p@sfile{#1}
			\def\@p@sbbfile{#1.bb}
			\@decmprtrue
		\fi}

\def\@p@@sfile#1{\@p@@sfigure{#1}}

\def\@p@@sbbllx#1{
		%\ps@typeout{bbllx is #1}
		\@bbllxtrue
		\dimen100=#1
		\edef\@p@sbbllx{\number\dimen100}
}
\def\@p@@sbblly#1{
		%\ps@typeout{bblly is #1}
		\@bbllytrue
		\dimen100=#1
		\edef\@p@sbblly{\number\dimen100}
}
\def\@p@@sbburx#1{
		%\ps@typeout{bburx is #1}
		\@bburxtrue
		\dimen100=#1
		\edef\@p@sbburx{\number\dimen100}
}
\def\@p@@sbbury#1{
		%\ps@typeout{bbury is #1}
		\@bburytrue
		\dimen100=#1
		\edef\@p@sbbury{\number\dimen100}
}
\def\@p@@sheight#1{
		\@heighttrue
		\dimen100=#1
   		\edef\@p@sheight{\number\dimen100}
		%\ps@typeout{Height is \@p@sheight}
}
\def\@p@@swidth#1{
		%\ps@typeout{Width is #1}
		\@widthtrue
		\dimen100=#1
		\edef\@p@swidth{\number\dimen100}
}
\def\@p@@srheight#1{
		%\ps@typeout{Reserved height is #1}
		\@rheighttrue
		\dimen100=#1
		\edef\@p@srheight{\number\dimen100}
}
\def\@p@@srwidth#1{
		%\ps@typeout{Reserved width is #1}
		\@rwidthtrue
		\dimen100=#1
		\edef\@p@srwidth{\number\dimen100}
}
\def\@p@@sangle#1{
		%\ps@typeout{Rotation is #1}
		\@angletrue
%		\dimen100=#1
		\edef\@p@sangle{#1} %\number\dimen100}
}
\def\@p@@ssilent#1{ 
		\@verbosefalse
}
\def\@p@@sprolog#1{\@prologfiletrue\def\@prologfileval{#1}}
\def\@p@@spostlog#1{\@postlogfiletrue\def\@postlogfileval{#1}}
\def\@cs@name#1{\csname #1\endcsname}
\def\@setparms#1=#2,{\@cs@name{@p@@s#1}{#2}}
%
% initialize the defaults (size the size of the figure)
%
\def\ps@init@parms{
		\@bbllxfalse \@bbllyfalse
		\@bburxfalse \@bburyfalse
		\@heightfalse \@widthfalse
		\@rheightfalse \@rwidthfalse
		\def\@p@sbbllx{}\def\@p@sbblly{}
		\def\@p@sbburx{}\def\@p@sbbury{}
		\def\@p@sheight{}\def\@p@swidth{}
		\def\@p@srheight{}\def\@p@srwidth{}
		\def\@p@sangle{0}
		\def\@p@sfile{} \def\@p@sbbfile{}
		\def\@p@scost{10}
		\def\@sc{}
		\@prologfilefalse
		\@postlogfilefalse
		\@clipfalse
		\if@noisy
			\@verbosetrue
		\else
			\@verbosefalse
		\fi
}
%
% Go through the options setting things up.
%
\def\parse@ps@parms#1{
	 	\@psdo\@psfiga:=#1\do
		   {\expandafter\@setparms\@psfiga,}}
%
% Compute bb height and width
%
\newif\ifno@bb
\def\bb@missing{
	\if@verbose{
		\ps@typeout{psfig: searching \@p@sbbfile \space  for bounding box}
	}\fi
	\no@bbtrue
	\epsf@getbb{\@p@sbbfile}
        \ifno@bb \else \bb@cull\epsf@llx\epsf@lly\epsf@urx\epsf@ury\fi
}	
\def\bb@cull#1#2#3#4{
	\dimen100=#1 bp\edef\@p@sbbllx{\number\dimen100}
	\dimen100=#2 bp\edef\@p@sbblly{\number\dimen100}
	\dimen100=#3 bp\edef\@p@sbburx{\number\dimen100}
	\dimen100=#4 bp\edef\@p@sbbury{\number\dimen100}
	\no@bbfalse
}
% rotate point (#1,#2) about (0,0).
% The sine and cosine of the angle are already stored in \sine and
% \cosine.  The result is placed in (\p@intvaluex, \p@intvaluey).
\newdimen\p@intvaluex
\newdimen\p@intvaluey
\def\rotate@#1#2{{\dimen0=#1 sp\dimen1=#2 sp
%            	calculate x' = x \cos\theta - y \sin\theta
		  \global\p@intvaluex=\cosine\dimen0
		  \dimen3=\sine\dimen1
		  \global\advance\p@intvaluex by -\dimen3
% 		calculate y' = x \sin\theta + y \cos\theta
		  \global\p@intvaluey=\sine\dimen0
		  \dimen3=\cosine\dimen1
		  \global\advance\p@intvaluey by \dimen3
		  }}
\def\compute@bb{
		\no@bbfalse
		\if@bbllx \else \no@bbtrue \fi
		\if@bblly \else \no@bbtrue \fi
		\if@bburx \else \no@bbtrue \fi
		\if@bbury \else \no@bbtrue \fi
		\ifno@bb \bb@missing \fi
		\ifno@bb \ps@typeout{FATAL ERROR: no bb supplied or found}
			\no-bb-error
		\fi
		%
%\ps@typeout{BB: \@p@sbbllx, \@p@sbblly, \@p@sbburx, \@p@sbbury} 
%
% store height/width of original (unrotated) bounding box
		\count203=\@p@sbburx
		\count204=\@p@sbbury
		\advance\count203 by -\@p@sbbllx
		\advance\count204 by -\@p@sbblly
		\edef\ps@bbw{\number\count203}
		\edef\ps@bbh{\number\count204}
		%\ps@typeout{ psbbh = \ps@bbh, psbbw = \ps@bbw }
		\if@angle 
			\Sine{\@p@sangle}\Cosine{\@p@sangle}
	        	{\dimen100=\maxdimen\xdef\r@p@sbbllx{\number\dimen100}
					    \xdef\r@p@sbblly{\number\dimen100}
			                    \xdef\r@p@sbburx{-\number\dimen100}
					    \xdef\r@p@sbbury{-\number\dimen100}}
%
% Need to rotate all four points and take the X-Y extremes of the new
% points as the new bounding box.
                        \def\minmaxtest{
			   \ifnum\number\p@intvaluex<\r@p@sbbllx
			      \xdef\r@p@sbbllx{\number\p@intvaluex}\fi
			   \ifnum\number\p@intvaluex>\r@p@sbburx
			      \xdef\r@p@sbburx{\number\p@intvaluex}\fi
			   \ifnum\number\p@intvaluey<\r@p@sbblly
			      \xdef\r@p@sbblly{\number\p@intvaluey}\fi
			   \ifnum\number\p@intvaluey>\r@p@sbbury
			      \xdef\r@p@sbbury{\number\p@intvaluey}\fi
			   }
%			lower left
			\rotate@{\@p@sbbllx}{\@p@sbblly}
			\minmaxtest
%			upper left
			\rotate@{\@p@sbbllx}{\@p@sbbury}
			\minmaxtest
%			lower right
			\rotate@{\@p@sbburx}{\@p@sbblly}
			\minmaxtest
%			upper right
			\rotate@{\@p@sbburx}{\@p@sbbury}
			\minmaxtest
			\edef\@p@sbbllx{\r@p@sbbllx}\edef\@p@sbblly{\r@p@sbblly}
			\edef\@p@sbburx{\r@p@sbburx}\edef\@p@sbbury{\r@p@sbbury}
%\ps@typeout{rotated BB: \r@p@sbbllx, \r@p@sbblly, \r@p@sbburx, \r@p@sbbury}
		\fi
		\count203=\@p@sbburx
		\count204=\@p@sbbury
		\advance\count203 by -\@p@sbbllx
		\advance\count204 by -\@p@sbblly
		\edef\@bbw{\number\count203}
		\edef\@bbh{\number\count204}
		%\ps@typeout{ bbh = \@bbh, bbw = \@bbw }
}
%
% \in@hundreds performs #1 * (#2 / #3) correct to the hundreds,
%	then leaves the result in @result
%
\def\in@hundreds#1#2#3{\count240=#2 \count241=#3
		     \count100=\count240	% 100 is first digit #2/#3
		     \divide\count100 by \count241
		     \count101=\count100
		     \multiply\count101 by \count241
		     \advance\count240 by -\count101
		     \multiply\count240 by 10
		     \count101=\count240	%101 is second digit of #2/#3
		     \divide\count101 by \count241
		     \count102=\count101
		     \multiply\count102 by \count241
		     \advance\count240 by -\count102
		     \multiply\count240 by 10
		     \count102=\count240	% 102 is the third digit
		     \divide\count102 by \count241
		     \count200=#1\count205=0
		     \count201=\count200
			\multiply\count201 by \count100
		 	\advance\count205 by \count201
		     \count201=\count200
			\divide\count201 by 10
			\multiply\count201 by \count101
			\advance\count205 by \count201
			%
		     \count201=\count200
			\divide\count201 by 100
			\multiply\count201 by \count102
			\advance\count205 by \count201
			%
		     \edef\@result{\number\count205}
}
\def\compute@wfromh{
		% computing : width = height * (bbw / bbh)
		\in@hundreds{\@p@sheight}{\@bbw}{\@bbh}
		%\ps@typeout{ \@p@sheight * \@bbw / \@bbh, = \@result }
		\edef\@p@swidth{\@result}
		%\ps@typeout{w from h: width is \@p@swidth}
}
\def\compute@hfromw{
		% computing : height = width * (bbh / bbw)
	        \in@hundreds{\@p@swidth}{\@bbh}{\@bbw}
		%\ps@typeout{ \@p@swidth * \@bbh / \@bbw = \@result }
		\edef\@p@sheight{\@result}
		%\ps@typeout{h from w : height is \@p@sheight}
}
\def\compute@handw{
		\if@height 
			\if@width
			\else
				\compute@wfromh
			\fi
		\else 
			\if@width
				\compute@hfromw
			\else
				\edef\@p@sheight{\@bbh}
				\edef\@p@swidth{\@bbw}
			\fi
		\fi
}
\def\compute@resv{
		\if@rheight \else \edef\@p@srheight{\@p@sheight} \fi
		\if@rwidth \else \edef\@p@srwidth{\@p@swidth} \fi
		%\ps@typeout{rheight = \@p@srheight, rwidth = \@p@srwidth}
}
%		
% Compute any missing values
\def\compute@sizes{
	\compute@bb
	\if@scalefirst\if@angle
% at this point the bounding box has been adjsuted correctly for
% rotation.  PSFIG does all of its scaling using \@bbh and \@bbw.  If
% a width= or height= was specified along with \psscalefirst, then the
% width=/height= value needs to be adjusted to match the new (rotated)
% bounding box size (specifed in \@bbw and \@bbh).
%    \ps@bbw       width=
%    -------  =  ---------- 
%    \@bbw       new width=
% so `new width=' = (width= * \@bbw) / \ps@bbw; where \ps@bbw is the
% width of the original (unrotated) bounding box.
	\if@width
	   \in@hundreds{\@p@swidth}{\@bbw}{\ps@bbw}
	   \edef\@p@swidth{\@result}
	\fi
	\if@height
	   \in@hundreds{\@p@sheight}{\@bbh}{\ps@bbh}
	   \edef\@p@sheight{\@result}
	\fi
	\fi\fi
	\compute@handw
	\compute@resv}

%
% \psfig
% usage : \psfig{file=, height=, width=, bbllx=, bblly=, bburx=, bbury=,
%			rheight=, rwidth=, clip=}
%
% "clip=" is a switch and takes no value, but the `=' must be present.
\def\psfig#1{\vbox {
	% do a zero width hard space so that a single
	% \psfig in a centering enviornment will behave nicely
	%{\setbox0=\hbox{\ }\ \hskip-\wd0}
	%
	\ps@init@parms
	\parse@ps@parms{#1}
	\compute@sizes
	%
	\ifnum\@p@scost<\@psdraft{
		%
		\special{ps::[begin] 	\@p@swidth \space \@p@sheight \space
				\@p@sbbllx \space \@p@sbblly \space
				\@p@sbburx \space \@p@sbbury \space
				startTexFig \space }
		\if@angle
			\special {ps:: \@p@sangle \space rotate \space} 
		\fi
		\if@clip{
			\if@verbose{
				\ps@typeout{(clip)}
			}\fi
			\special{ps:: doclip \space }
		}\fi
		\if@prologfile
		    \special{ps: plotfile \@prologfileval \space } \fi
		\if@decmpr{
			\if@verbose{
				\ps@typeout{psfig: including \@p@sfile.Z \space }
			}\fi
			\special{ps: plotfile "`zcat \@p@sfile.Z" \space }
		}\else{
			\if@verbose{
				\ps@typeout{psfig: including \@p@sfile \space }
			}\fi
			\special{ps: plotfile \@p@sfile \space }
		}\fi
		\if@postlogfile
		    \special{ps: plotfile \@postlogfileval \space } \fi
		\special{ps::[end] endTexFig \space }
		% Create the vbox to reserve the space for the figure.
		\vbox to \@p@srheight sp{
		% 1/92 TJD Changed from "true sp" to "sp" for magnification.
			\hbox to \@p@srwidth sp{
				\hss
			}
		\vss
		}
	}\else{
		% draft figure, just reserve the space and print the
		% path name.
		\if@draftbox{		
			% Verbose draft: print file name in box
			\hbox{\frame{\vbox to \@p@srheight sp{
			\vss
			\hbox to \@p@srwidth sp{ \hss \@p@sfile \hss }
			\vss
			}}}
		}\else{
			% Non-verbose draft
			\vbox to \@p@srheight sp{
			\vss
			\hbox to \@p@srwidth sp{\hss}
			\vss
			}
		}\fi	



	}\fi
}}
\psfigRestoreAt
\let\@=\LaTeXAtSign

astro-ph0402557.extracted_bib
\begin{thebibliography}{}
\def\ref#1\par{\parshape=2 0in 14.5cm 1cm 13.5cm {#1} \par}
\parskip=0pt
\parindent=0pt

\bibitem[]{}
Chakrabarti S.K. 1989, ApJ, 347, 365

\bibitem[]{}
Chakrabarti S.K. 1996a, ApJ, 464, 664

\bibitem[]{}
Chakrabarti S.K. 1996b, MNRAS, 283, 325

\bibitem[]{}
Chakrabarti S.K. 1999, A\&A 351, 185

\bibitem[]{}
Chakrabarti S.K., Manickam, S.G. 2000, ApJ 531, L41 (CM00)

\bibitem[]{}
Chakrabarti S.K., Molteni, D.G. 1995, MNRAS, 272, 80

\bibitem[]{} 
Chakrabarti, S.K.,  Nandi, A., Manickam, S. G.,  Mandal, S.  Rao, A.R. 2002, ApJ 579, L21

\bibitem[]{}
Chanmugam, G., Langer, S.H., Shaviv, G. 1985, ApJ, 288, 428

\bibitem[]{}
Dermer, C.D., Liang, E.P., Canfield, E. 1991, ApJ, 369, 410

\bibitem[]{}
Dhawan, V., Mirabel, I. F., Rodriguez, L. F. 2000, ApJ 543 373

\bibitem[]{}
Halpern, J.P., Marshall, H.L. 1996,  ApJ, 464, 760

\bibitem[]{}
Landau, L.D., Lifshitz, E.D. 1959, Fluid Mechanics (New York:Pergamon)

\bibitem[]{}
Langer, S.H., Chanmugam, G., Shaviv, G. 1981, ApJ, 245, L23

\bibitem[]{} 
Molteni, D., Acharyya, K., Kuznetsov, O., Bisikalo, D.,  Chakrabarti, S.\ K. 2001, ApJL, 563, L57

\bibitem[]{} 
Molteni, D., Sponholz, H.,  Chakrabarti, S.\ K. 1996, ApJ, 457, 805  (MSC96)

\bibitem[]{} 
Molteni, D., Ryu, D., Chakrabarti, S.\ K. 1996, ApJ, 470, 460

\bibitem[]{}
Paczy\'nski, B., Wiita, P.J., 1980, A\&A, 88, 23 

\bibitem[]{}
Rao, A.R., Naik, S., Vadawale, S.V., Chakrabarti, S.K. 2000, A\&A 360, L25

\bibitem[]{}
Rees, M.J. 1984, Ann. Rev. Astron. Astrophys., 22, 471

\bibitem[]{}
Rybicki, G., Lightman, A.P., 1979, Radiative Processes in Astrophysics (New York:Wiley Interscience)

\bibitem[]{}
Ryu, D., Chakrabarti, S.K., Molteni, D. 1997, ApJ, 474, 378 (RCM97)

\bibitem[]{}
Shaham, J. 1987 in The Origin and Evolution of Neutron Stars, D.J. Helfand and J.-H. Huang (eds.), 347

\bibitem[]{}
Smith, D. M., Heindl, W. A., Swank, J. H., 2002, ApJ, 569, 362

\end{thebibliography}