\documentstyle[epsf]{l-aa}
%\documentstyle[referee]{l-aa}

\hyphenation{pass-band pass-bands Eq two scale-lengths gal-axies}
\def\ha{\relax \ifmmode {\rm H}\alpha\else H$\alpha$\fi}
\def\hi{\relax \ifmmode {\rm H\,{\sc i}}\else H\,{\sc i}\fi}
\def\hii{\relax \ifmmode {\rm H\,{\sc ii}}\else H\,{\sc ii}\fi}
\def\rq{$r^{1/4} $ }
\def\rh{$r^{1/2} $ }
\def\muo{\relax \ifmmode \mu_0\else $\mu_0$\fi}
\def\mue{\relax \ifmmode \mu_{\rm e}\else $\mu_{\rm e}$\fi}
\def\re{\relax \ifmmode r_{\rm e}\else $r_{\rm e}$\fi}
\def\magarc{mag arcsec$^{-2}$}
\def\boundboxo{\epsfbox[30 190 538 540]}
\def\boundboxt{\epsfbox[10 190 518 540]}

\def\abst{
 The stellar and dust content of spiral galaxies as function of radius
has been investigated using near-infrared and optical broadband surface
photometry of 86 face-on spiral galaxies.  Colors of galaxies correlate
with the azimuthally averaged local surface brightness both within and
among galaxies, with the lower surface brightness regions being bluer. 
The colors formed from different passband combinations correlate
strongly indicating that they probably arise from the same physical
process. 

A 3D radiative transfer model was developed to calculate the effect of
dust absorption and scattering on the luminosity and color profiles of
galaxies.  Stellar synthesis models were used to investigate the effects
of the star formation history and the metallicity on the broadband color
profiles.  Combining all optical and near-infrared data shows that the
color gradients in this sample of face-on galaxies are best explained by
a combined stellar age and metallicity gradient across the disk, with
the outer regions being on average younger and of lower metallicity. 
Dust reddening probably plays only a minor role, as the dust models
cannot produce reddening profiles that are compatible with the
observations. 

The observed color differences implicate substantial $M/L_\lambda$
differences, both within galaxies and among galaxies.  The variations
are such that the ``missing light'' problem derived from rotation
fitting becomes even worse. Late-type galaxies (T$\ge$6) have lower
metallicities and are often of younger average age than earlier types
and have therefore an entirely different $M/L_\lambda$ in most
passbands.
 %Late-type galaxies should be used with caution in the
%Tully-Fisher relation.
 The near-infrared passbands are recommended for studies where the
$M/L_\lambda$ ratios should not vary too much. 

\keywords{dust, extinction -- Galaxies: evolution --Galaxies:
fundamental parameters -- Galaxies: photometry -- Galaxies: spiral --
Galaxies: stellar content}
}


\begin{document}
\onecolumn

 \thesaurus{03(09.04.1; 11.05.2; 11.06.2; 11.16.1; 11.19.2; 11.19.5)}
 \title{
Near-infrared and optical broadband surface photometry of 86 face-on
disk dominated galaxies. 
 \thanks{Based on observations with the Jacobus Kapteyn Telescope and the
Isaac Newton Telescope operated by the Royal Greenwich Observatory at
the Observatorio del Roque de los Muchachos of the Instituto de
Astrof\'{\i}sica de Canarias with financial support from the PPARC (UK) and NWO
(NL) and with the UK Infrared Telescope at Mauna Kea operated by the
Royal Observatory Edinburgh with financial support of the PPARC.
\newline
The luminosity/color profiles of all galaxies are available in
electronic form at the CDS via anonymous ftp 130.79.128.5 as part of Paper~I.
 }
} 
 \subtitle{IV. Using color profiles to study stellar and dust content of galaxies.}
 \author{Roelof S.~de Jong \inst{ } }
 \offprints{
R.S.~de Jong, University of Durham, Dept.~of
Physics, South Road, Durham DH1 3LE, United Kingdom, e-mail:
R.S.deJong@Durham.ac.uk} 
 \institute{
Kapteyn Astronomical Institute, P.O.Box 800, NL-9700 AV
Groningen, The Netherlands}
 \date{received, accepted}
 \maketitle
 \markboth{R.S.~de Jong}{Near-IR and
 optical broadband observations of 86 spirals. IV Color profiles}



\begin{abstract}
\abst
\end{abstract}

\appendix
\section{Monte Carlo radiative transfer simulations of
light and dust in exponential disks}
\label{apmodel}

The modeling of dust extinction in extragalactic systems has a long
history. The effects of scattering are often ignored or are assumed to be
no more than a scaling factor (Disney et al.~\cite{DDP89};
Huizinga~\cite{HuiPhD}), which is correct as long as one is not looking
at wavelength dependent effects. Therefore, scattering is included in
most studies investigating the wavelength dependent effects of
extinction (e.g.\ Kylafis \& Bahcall~\cite{KylBah87}; Bruzual et
al.~\cite{Bru88}, Witt et al.~\cite{Witt92}, Byun et al.~\cite{Byun94})

The Monte Carlo approach followed here is a somewhat ``brute force'' method,
but it has the advantage that it can be used for luminous and dust
geometries containing little symmetry. Even though smooth light and dust
configurations are applied here, this method can easily be extended to
include dust clouds and spiral structure.

 \subsection{The mathematical method}
 A Monte Carlo radiative transfer code was used to calculate the light
and color distribution of disk galaxies as seen by a distant observer.
The Monte Carlo principle as applied to extinction in gaseous nebulae
is described in detail by Witt (\cite{Witt77}).  In the computer model
the paths of many ($\sim10^6$) photons were followed as they
traveled through the absorbing and scattering dusty medium.  At
great distance the photons were collected to produce a galaxy image
which was used for further study. 

The trajectory of a photon through a dusty medium is determined by
random processes.  One can characterize these processes by a probability
function $p(x)$ on an interval $(a,b)$ such that
 \begin{equation}
\int_{a}^{b} p(\xi ) d\xi = 1. 
 \end{equation}
Using a random number generator, which produces a random number
$R$ distributed uniformly in the interval $0 \!\le \!R \!\le \!1$, one can
simulate an event with frequency $p(x)dx$ in the interval $(x,x + dx)$
by requiring 
 \begin{equation} 
\int_{a}^{x} p(\xi ) d\xi = R.
\label{Rint}
 \end{equation}
 This equation was used to simulate the birthplace and direction of
photons, as well as the scattering properties of the dust. In the
following sections, each $R$ will denote a new random number in the interval
$0\!\le \!R \!\le \!1$.


\subsection{The creation of photons}

%\input{coorfig}
\begin{figure*}
 \mbox{\epsfxsize=17.5cm\epsfbox[10 340 500 630]{coorfig2.ps}}
 \caption[]{
 The definition of the coordinate system.  A trajectory with a
scattering is indicated for one step size.  A ``photon beam''
starts at $\overline{x}_{\rm old}$ in the direction defined by
($\omega$,$\theta$), scatters off a dust particle over an angle $\alpha$,
and ends at $\overline{x}_{\rm new}$. 
}
 \label{coorfig}
 \end{figure*}


This section describes the creation of photons and the random processes
involved in this creation. Rather than creating photons in space with
a certain density distribution such as in a real galaxy, the photons
were instead created uniformly in space and given an initial intensity
weight to produce the exponential light profile. These initial weights
were reduced by absorption as the photon moved through the dusty
medium.

The models consisted of three-dimensional distributions of stellar light
and dust that were chosen independently.  For the stellar light, an
axisymmetric disk-like distribution was used, with an exponential
intensity behavior in both radial ($r$) and vertical ($z$) directions:
 \begin{equation}
I(r,z,\phi) = {\rm e}^{-(r/h_{\rm s}+|z|/z_{\rm s})},
\label{inint}
 \end{equation}
 where $h_{\rm s}$ and $z_{\rm s}$ are the scalelength and the
scaleheight of the stellar distribution respectively, with the
coordinate system is as defined in Fig.~\ref{coorfig}.  The luminosity was
truncated at seven scalelengths and heights.  The exponential behavior
of the radial light distribution of disks in spiral galaxies is well
established, but the vertical light distribution is more controversial,
because extinction effects make measurements difficult.  I have used
the exponential law vertically instead of (for instance) the sech or
sech$^2$ laws (van der Kruit~\cite{Kru88}), because the nearly
unobscured light distribution as obtained with near-IR observations
often can be well described by such an exponential law (Wainscoat et
al.~\cite{Wai89}; Aoki et al.~\cite{Aoki91}).  It is not expected that
the results obtained here will change significantly if other
plausible vertical light distributions are used.  Each ``photon beam''
that was created received an initial weight according to the
$I(r,z,\phi)$ of Eq.~(\ref{inint}). 

Photons were created at a certain position in the model galaxy using the
distribution functions
 \begin{equation}
r = 7 R h_{\rm s} \ \ \ \ 
z = 7 (R-0.5)2 z_{\rm s} \ \ \ \ 
\phi = 2 R \pi,
\label{inpos}
 \end{equation}
and thus their initial position in cartesian coordinates was
 \begin{equation}
\overline{x} = \left( \begin{array}{c} r \cos(\phi) \\
                                       r \sin(\phi) \\
                                       z
                      \end{array} \right).
 \end{equation}
 The distribution functions of Eq.~(\ref{inpos}) create photons
uniformly in cylindrical coordinates, but not in cartesian space.  The
density of created photons is much higher near the center than in the
outer regions and a correction factor is needed.  To still get an
exponential behavior in flux, an additional weight factor of $2 \pi r$
was given to the intensity of each created ``photon beam''.  These
distribution functions were chosen, because the final model images were
azimuthally averaged just as the real observations and in the absence of
dust give these distribution functions an equal number of photons at
each radius in the face-on case. 


At creation, the initial flight direction of the photons was specified
by the functions
 \begin{equation}
\omega = (2 R - 1) \pi \ \ \ \  
\theta = (R - 0.5) \pi .
\label{indir}
 \end{equation}
The corresponding directional cosines are
 \begin{equation}
\Delta \overline{x} = \left( \begin{array}{c} \cos(\omega) \cos(\theta) \\
                                              \sin(\omega) \cos(\theta) \\
                                              \sin(\theta) 
                             \end{array} \right).
\label{dircos}
 \end{equation}
 Again, since this would not give a uniform distribution of flux
density in each direction, an extra weight of the form $\cos(\theta)$
was introduced. The initial intensity of a ``photon beam'' at creation
was thus
 \begin{equation}
 I_0(r,z,\phi,\omega,\theta) = 2 \pi r \cos(\theta) {\rm e}^{-(r/h_{\rm s}+|z|/z_{\rm s})}.
 \end{equation}

\subsection{The dust properties}

Let us first define the notation for the basic equations of radiative
transfer.  A beam of light with intensity $I$ is followed while
traveling along a straight line through a dusty medium.  Since this medium
removes a fraction $\kappa\, d s$ of the incident intensity for each
distance $d s$ traveled, the full equation is
 \begin{equation}
\frac{d I(s)}{d s}  =  -\kappa(s) I(s) + J(s),
 \end{equation}
 where $J(s)$ are the new sources of intensity at point $s$ in the
travel direction. $\kappa$ is the extinction coefficient and
consists of an absorption and a scattering part, and their relative
importance is normally expressed by the albedo ($a$)
 \begin{eqnarray}
\kappa & = & \kappa_{a}+\kappa_{s} \\
a & = & \kappa_{s}/\kappa.
 \end{eqnarray}
  While following the light beam through the dusty medium, no new photons
were added to the initial beam of photons that was created with
intensity $I(0)$, because they are accounted for by the repeating Monte
Carlo principle, and $J(s) \!= \!0$.  For this single ``photon beam''
the radiative transfer equation is now
 \begin{equation}
\frac{d I(s)}{d s}  =  -\kappa(s) I(s)
 \end{equation}
with formal solution
 \begin{equation}
I(x)  =  I(0) {\rm e}^{-\int_{0}^{x} \kappa(s) ds}.
\label{inteq}
 \end{equation}
The integral in Eq.~(\ref{inteq}) defines optical thickness
 \begin{equation}
\tau \equiv \int_{0}^{x} \kappa(s) \,ds.
 \end{equation}
 Since the beam was followed while it scattered through the medium (so
no longer on a straight line), there were no losses by scattering along
the path of the beam.  One still can use Eq.~(\ref{inteq}), but with
$\kappa_{a}(s) \!= \! (1 \!- \!a)\kappa(s)$ instead of $\kappa(s)$ as
long as the line integral is taken along the path traveled.  The
extinction coefficient and the albedo are wavelength dependent, which
creates reddening of stars and could cause the observed color gradients
in the galaxies. 

To calculate the absorption along the path traveled, the dust
distribution must be defined.  The dust extinction coefficient in the
models was taken to have a similar distribution as the stellar light,
but with a scalelength and scaleheight that were chosen independently:
 \begin{equation}
\kappa_{\lambda}(r,z,\phi) = \kappa_{0,\lambda} {\rm e}^{-(r/h_{\rm d}+|z|/z_{\rm d})}
 \end{equation}
 where $\kappa_{\lambda}(r,z,\phi)$ is the local extinction coefficient at
wavelength $\lambda$ and $\kappa_{0,\lambda}$ is the extinction
coefficient at $(r\!=\!0, z\!=\!0)$.  In 
this article {\em the} optical depth of a system is denoted by the
integration of $\kappa_V$ ($\kappa$ in the $V$-passband) along the
symmetry axis from $(r\!=\!0,
z\!=\!-\infty)$ to $(r\!=\!0, z\!=\!\infty)$
 \begin{equation}
\tau_{0,V} = \int_{-\infty}^{\ \infty} \kappa_{0,V} {\rm e}^{-|z|/z_{\rm d}} \,dz 
      = 2 \kappa_{0,V} z_{\rm d}.
 \label{optdepth}
 \end{equation}
 Using Eq.~(\ref{inteq}) one can calculate that a point source located
behind the center of the galaxy will have suffered an extinction in the
$V$-passband for an observer located at its other pole of
 \begin{equation}
I(z=\infty)  =  I(0) {\rm e}^{-\tau_{0,V}},
\label{exflux}
 \end{equation}
assuming that all absorbed {\em and} scattered photons are lost for a
point source.

{\tabcolsep=4mm
\begin{table}
\begin{center}
 \begin{tabular}{cccc}
\hline
\hline
\ \ passband & $\tau_\lambda/\tau_V$ & $a_\lambda$ & $g_\lambda$\\
\hline
$U$     & 1.531 & 0.68  & 0.67\\
$B$	& 1.324	& 0.66	& 0.59\\
$V$	& 1.000 & 0.60	& 0.50\\
$R$	& 0.748 & 0.53	& 0.40\\
$I$	& 0.482 & 0.45	& 0.29\\
$H$	& 0.175 & 0.28	& 0.04\\
$K$	& 0.112	& 0.20	& 0.00\\
\hline
\hline
\end{tabular}
 \caption[]{Values used in the scattering model for the dust properties
relative extinction ($\tau_\lambda/\tau_V$), albedo ($a_\lambda$) and
scattering asymmetry ($g_\lambda$) as function of photometric passband. }
\label{dusprop}
\end{center}
\end{table}
}


Three parameters are needed to describe dust properties,
the relative extinction $(\tau_\lambda/\tau_V)$, the albedo
($a_\lambda$) and the scattering phase function $\Phi_\lambda(\alpha)$,
and all are dependent on wavelength $(\lambda)$.  As there are almost no direct
measurements of the dust properties in other galaxies, I decided to use
the (also poorly determined) Galactic values.  Studies by Knapen et
al.\ (\cite{Kna91}) and Jansen et al.~(\cite{Jan94}) seem to indicate
that at least the Galactic extinction curve is applicable to some other
galaxies, the other two dust properties have never been measured in
extragalactic systems.

For relative extinction, the Galactic extinction properties of Rieke
\& Lebofsky (\cite{RieLeb85}) were used. The values for the albedo
were drawn from Bruzual et al.~(\cite{Bru88}), as were the values for
the scattering asymmetry parameter $g_\lambda$. This asymmetry
parameter enters in the scattering phase function suggested by Henyey
\& Greenstein (\cite{Henyey})
 \begin{eqnarray}
\lefteqn{\Phi_\lambda(\cos(\alpha),g_\lambda) = }  \nonumber\\
& \ \ & [(1- g_\lambda^2)/4\pi] (1+g_\lambda^2-2g_\lambda\cos(\alpha))^{-3/2}
\label{HenGre}
\end{eqnarray}
 where $\alpha$ is the scattering angle between the incident and the
deflected photon.  The function is such that $g_\lambda \!= \!\langle
\cos(\alpha) \rangle$, and thus $-1 \!\leq \!g_\lambda \!\leq \!1$. 
Forward scattering dominates for $g_\lambda \!> \!0$ and $g_\lambda \!=
\!0$ results in isotropic scattering.
 %Note that Eq.~(\ref{HenGre}) is in fact a distribution function of
%$\cos(\alpha)$ and not of $\alpha$ as suggested in the original paper
%by Henyey \& Greenstein (\cite{Henyey}). 
 To determine a random deflection angle
for a photon with the probability function $\Phi_\lambda$ 
Eq.~(\ref{Rint}) can be used and one finds (Witt \cite{Witt77})
 \begin{eqnarray}
\alpha & = & \arccos(\{(1+g_\lambda^2)- \nonumber \\
 & & [(1-g_\lambda^2)/(1-g_\lambda+2g_\lambda R)]^2\}/2g_\lambda).
\label{HenGreWit}
 \end{eqnarray} 
 The values used for $\tau_\lambda/\tau_V$, $a_\lambda$, and $g_\lambda$
for the different photometric passbands are listed in
Table~\ref{dusprop}.

%\onecolumn

\subsection{The numerical method}

Each photon beam that was created using Eqs.~(\ref{inpos}) and
(\ref{indir}) was stepped through the dusty medium. The step size was
determined by the local dust density and was chosen in such
way that it would give an optical thickness (absorption plus
scattering) of 0.03 if there was a constant dust density along the
step. This means that the chance of multiple scattering is of order
$(0.03 a_\lambda)^2$ along a step. This is negligible, as it should be,
because only single scattering was incorporated along each step.  It is
trivial to show that the step size has to be of order
 \begin{equation}
\Delta s = 0.03 / \kappa(\overline{x})
         = 0.03 / (\kappa_0 {\rm e}^{-(r/h_{\rm s}+|z|/z_{\rm s})})
\label{step}
 \end{equation}
and the new position is (if there is no scattering)
 \begin{equation}
\overline{x}_{\rm new} = \overline{x}_{\rm old} + \Delta \overline{x} \Delta s.
 \end{equation}

For the actual absorption along a step, one has to use Eq.~(\ref{inteq})
with $\kappa_a \!= \!\kappa(1-a)$ instead of $\kappa$.  The integral was
at each step approximated by four-point Gauss-Legendre quadrature:
% \begin{eqnarray}
%\tau_{\Delta s} & =       & \int_{\overline{x}}^{\overline{x}+\Delta \overline{x} \Delta s} \kappa(\overline{s})d\overline{s} \\
%                & \approx &\Delta s [\kappa(\overline{x})
%                            + 4 \kappa(\overline{x} +\Delta \overline{x} \Delta s/2) 
%                            +\kappa(\overline{x} +\Delta \overline{x} \Delta s) ]/6. \nonumber
% \end{eqnarray}
% \begin{eqnarray}
%\tau_{\Delta s} & =       & \int_{\overline{x}}^{\overline{x}+\Delta \overline{x} \Delta s} \kappa_a(\overline{s})\,d\overline{s} \\
%                & \approx & 0.5 \Delta s \sum_{j=1}^4 w_j
%\,\kappa(\overline{x}+0.5(1+x_j)\Delta\overline{x})(1-a)
%\label{GauLeg}
% \end{eqnarray}
 \begin{equation}
\tau_{\Delta s}  =      \int_{\overline{x}}^{\overline{x}+\Delta \overline{x} \Delta s} \kappa_a(\overline{s})\,d\overline{s} 
                \ \ \approx \ \ 0.5 \Delta s \sum_{j=1}^4 w_j
\,\kappa(\overline{x}+0.5(1+x_j)\Delta\overline{x})(1-a),
\label{GauLeg}
 \end{equation}
where the values for $w_j$ and $x_j$ can be found in any handbook on
numerical analysis.
%with
% \begin{eqnarray}
% w_{1,2} = 0.3478548451 &{\rm and}& \overline{x}_{1,2} = \overline{x}+\Delta\overline{x} \Delta s (1\pm0.8611363116)/2  \\
% w_{3,4} = 0.6521451549 &{\rm and}& \overline{x}_{3,4} = \overline{x}+\Delta\overline{x} \Delta s (1\pm0.3399810436)/2 
% \end{eqnarray}
 When there was no scattering during the step, the
$\kappa(\overline{x}_4)$ calculated in Eq.~(\ref{GauLeg}) was used to
determine the step size with Eq.~(\ref{step}) in the next step.

%The weight of a photon was reduced during a step by a factor
%e$^{-\tau}$, i.e. $w(\overline{x}_{\rm new}) = w(\overline{x}_{\rm old)
%e$^{-\tau}$.

If scattering did occur, the step proceeded slightly differently. The
chance of scattering during a step can be simulated by requiring
 \begin{equation}
{\rm e}^{-\tau_{\Delta s} a_\lambda} < R .
\label{scatchan}
 \end{equation}
In that case the directional cosines of Eq.~(\ref{dircos}) were changed
by an angle given by Eq.~(\ref{HenGreWit}). This left another angle of
freedom for the azimuthal change of the direction which was generated by
 \begin{equation}
\beta = 2\pi R.
 \end{equation}
Rotating the initial directional angles $\omega,\theta$ over the
scattering angles $\alpha,\beta$ gives the new directional cosines
 \begin{equation}
\Delta \overline{x}_{\rm new} =\left( \begin{array}{rr}
\cos(\omega) \cos(\theta) \cos(\beta)&+( \sin(\omega) \cos(\alpha)+\cos(\omega) \sin(\theta) \sin(\alpha) ) \sin(\beta) \\
\sin(\omega) \cos(\theta) \cos(\beta)&+(-\cos(\omega) \cos(\alpha)+\sin(\omega) \sin(\theta) \sin(\alpha) ) \sin(\beta) \\
\sin(\theta)              \cos(\beta)&                                          -    \cos(\theta) \sin(\alpha) \sin(\beta) \\
                     \end{array} \right).
 \end{equation} 
 % \begin{equation} 
%\Delta \overline{x}_{\rm new} =\left( \begin{array}{rrr} 
%\cos(\omega) \cos(\theta)& \sin(\omega) &\cos(\omega) \sin(\theta) \\ 
%\sin(\omega) \cos(\theta)&-\cos(\omega) &\sin(\omega) \sin(\theta) \\ 
%\sin(\theta)             &             &            -\cos(\theta) \\ 
%                     \end{array} \right)
%\left( \begin{array}{r} 
%\cos(\beta)\\ 
%\cos(\alpha)\sin(\beta)\\
%\sin(\alpha)\sin(\beta)\\ 
%\end{array} \right) 
% \end{equation} 
%\begin{equation} 
%\Delta \overline{x}_{\rm new} =\left(\begin{array}{rr} 
%c(\omega) c(\theta) c(\beta)&+( s(\omega) c(\alpha)+c(\omega) s(\theta) s(\alpha) ) s(\beta) \\ 
%s(\omega) c(\theta) c(\beta)&+(-c(\omega) c(\alpha)+s(\omega) s(\theta) s(\alpha)) s(\beta) \\ 
%s(\theta)              c(\beta)&                                         -    c(\theta) s(\alpha) s(\beta) \\ 
%                    \end{array} \right) 
 % \end{equation}
 When Eq.~(\ref{scatchan}) was satisfied and scattering occurred in a step,
the photon proceeded a fraction $f \!= \!R$ of step size $\Delta s$ in the
old direction, before following the new direction
 \begin{equation}
\overline{x}_{\rm new} = \overline{x}_{\rm old} + (f \Delta
\overline{x}_{\rm old}  + (1-f) \Delta \overline{x}_{\rm new}) \Delta s.
 \end{equation}

\subsection{Projection on the sky}

The steps in the previous paragraph were repeated until
$|\overline{x}_{\rm new}| \!> \!10 h_{\rm s}$. All $\tau_{\Delta s}$ of
the different steps were added, and the final intensity of the photon
beam was $I_{\rm end} \!= \!I_0 {\rm e}^{-\tau_{\rm tot}}$. The final
intensities were projected on the sky as if the galaxy was being observed
from infinity. Due to the axisymmetric nature of the models, one can
ignore the $\omega$ dependence of the exit direction and rotate all
photons over angle $\omega$ as if they leave the galaxy in the same
direction. The projection on the $(y_1,y_2)$-plane becomes
 \begin{equation}
\overline{y} = \left( \begin{array}{rrr}
-\sin(\omega)            & \cos(\omega)            & 0\\
-\sin(\theta)\cos(\omega)&-\sin(\theta)\sin(\omega)&\cos(\theta)
\end{array} \right) \overline{x}
 \end{equation}
 This procedure of creating photons, stepping through the medium and
projecting the exiting photons on the sky was repeated for two million
photons.  The photons were binned into pixel images in $(y_1,y_2)$
direction and the binning in viewing angle $\theta$ direction resulted
in ten model images from edge-on to face-on.  The binning of photons in
the $\theta$ direction was in equal steps of $b/a$, where $b/a$ is the
axial ratio of an inclined circle.  This has the advantage that each
image has approximately equal flux, at least in the case of no
extinction.  For the analyses discussed in this paper only the face-on
images were used, which contained the photons that exited with an angle
between 90\degr\ and 71.8\degr\ (i.e.\ $1 \!\ge \!b/a \!> \!0.95$). 

%\twocolumn
\subsection{Testing and the results}

The computer program was tested extensively to ensure that the results
were reliable.  First, a point source emitting all its photons to a pole
was placed at the center of the galaxy and the exiting flux was
correctly described by Eq.~(\ref{exflux}).  Tests showed that in the
case of no absorption all images at different viewing angle had an equal
flux.  Tests were made to ensure that the scattering angle phase
function was produced correctly
 %(Fig.~\ref{phidis}) 
 and that the angle between the incident and the deflected photon was
indeed the required angle.  One of the most important tests was the
independence of the number of scatterings on step size.  If the step
size had been too large or if the calculation $\tau_{\Delta s}$ had been
incorrect this could not have been the case.  The model profiles were
tested against the analytic profiles of DDP for the cases with and
without scattering.  In the case of no scattering (all scattered photons
were removed to mimic absorption) the agreement was perfect. 
Figure~\ref{duscolcol} shows an excellent agreement between the two
models even with scattering, as long as the extinction is high.  All
test results were correct to within the statistical noise. 

For the extraction of the (color)profiles, the same programs as in
Paper~I were used and the model results can be compared directly with
the data. Once a $\tau_{0,V}$ and a set of scalelengths and heights have
been chosen for a particular model, $\kappa_{0,\lambda}$ is determined
by Eq.(~\ref{optdepth}) and the corresponding optical depths for the other
passbands are determined by Table~\ref{dusprop}.  This model is also
useful to study inclination-dependent effects of extinction. 

\begin{figure*}[tbh]
 \mbox{\epsfxsize=8.8cm\epsfbox[20 122 575 670]{bvrk1.ps}}
 \mbox{\epsfxsize=8.8cm\epsfbox[20 122 575 670]{bivk1.ps}}
 \mbox{\epsfxsize=8.8cm\epsfbox[20 163 575 670]{bvrk2.ps}}
\ \ \ 
 \mbox{\epsfxsize=8.8cm\epsfbox[20 163 575 670]{bivk2.ps}}
 \caption[]{
 Color--color plots resulting from dust models of face-on galaxies.  The
color zero-points are arbitrary (depending on the color of the
underlying stellar population) and are indicated by the filled circles. 
The centers of the model galaxies are indicated by the open circles and the
colors are followed in radial direction for four and a half scale
lengths.  The solid lines are the results from the Monte Carlo
simulations, described in Appendix~\ref{apmodel}.  The dotted lines are
the results from the DDP models, using the Galactic absorption law
rather than the extinction law between the different passbands.  The
reddening vectors in the top right corners indicate the Galactic
extinction law (i.e.\ Screen model). 
\vspace{-2ex}
 }
 \label{duscolcol}
\end{figure*}



\begin{thebibliography}{}
%
%\bibitem[1978]{Aar78}
%Aaronson M. 1978, ApJ 221, L103
%
%\bibitem[1992]{AndKru92}
%Andredakis Y.C., van der Kruit P.C. 1992, A\&A 265, 396
%
\bibitem[1991]{Aoki91}
Aoki T.E., Hiromoto N., Takami H., Okamura S. 1991, PASJ 43, 755
%
%\bibitem[1944]{Baa44}
%Baade W. 1944, ApJ 100, 137
%
%\bibitem[1987]{BegPhD}
%Begeman K. 1987, Ph.D.~Thesis, University of Groningen, The Netherlands
%
%\bibitem[1993]{Ber93}
%Bershady M.A. 1993, PASP 105, 1028
%
%\bibitem[1979]{Bes79} %Don't delete!
%Bessell M.S. 1979, PASP 91, 589
%
%\bibitem[1994]{Blo94}
%Block D.L., Witt A.N., Grosb\o l P., Stockton A., Moneti A. 1994, A\&A
%288, 383
%
%\bibitem[1978]{BosPhD}
%Bosma A. 1978, Ph.D.~Thesis, University of Groningen, The Netherlands
%
%\bibitem[1993]{Bre93}
%Bressan A., Fagotto F., Bertelli G., Chiosi C. 1993, A\&AS 100, 647
%
%\bibitem[1992]{BroPhD}
%Broeils A. 1992, Ph.D.~Thesis, University of Groningen, The Netherlands
%
%\bibitem[1993]{BruCha93}
%Bruzual G.A., Charlot S. 1993, ApJ 405, 538
%
%\bibitem[1996]{BruCha96}
%Bruzual G.A., Charlot S. 1996, ApJ in preparation ({\bf BC96})
%
\bibitem[1988]{Bru88}
Bruzual G.A., Magris G.C., Calvet N. 1988, ApJ 333, 673
%
%\bibitem[1984]{BurHei84}
%Burstein D., Heiles C. 1984, ApJS 54, 33
%
%%\bibitem[1992]{BusSta92}
%%Bushouse H.A., Stanford S.A. 1992, ApJS 79, 213
%
%\bibitem[1989]{Buz89}
%Buzoni A. 1989, ApJS 71, 817
%
\bibitem[1994]{Byun94}
Byun Y.I., Freeman K.C., Kylafis N.D. 1994, ApJ 432, 114
%
%\bibitem[1991]{ChaBru91}
%Charlot S., Bruzual G.A. 1991, ApJ 367, 126
%
%\bibitem[1995]{CWB95}
%Charlot S., Worthey G., Bressan A. 1995, ApJ submitted
%
%\bibitem[1987]{Cor87}
%Cornell M., E., Aaronson M., Bothun G., Mould J. 1987, ApJS 64, 507
%
%\bibitem[1995]{deJPhD}
%de Jong R.S. 1995, Ph.D.~Thesis, University of Groningen, The Netherlands
%
%\bibitem[1996a]{deJ2}
%de Jong R.S. 1996a, A\&A in press ({\bf Paper~II})
%
%\bibitem[1996b]{deJ3}
%de Jong R.S. 1996b, A\&A in press ({\bf Paper~III})
%
%\bibitem[1996c]{deJeso} de Jong R.S.  1996c, In: Spiral Galaxies in
%the near-infrared, eds.\ D.~Minniti, H.-W.~Rix (ESO Workshop series), in
%press
%
%\bibitem[1994]{deJ1}
%de Jong R.S., van der Kruit P.C. 1994, A\&AS 106, 451 ({\bf Paper~I})
%
%\bibitem[1991]{rc3}
%de Vaucouleurs G., de Vaucouleurs A., Corwin H.G., Buta R.J. et al. 1991,
%Third Reference Catalog of Bright Galaxies (Springer-Verlag, New York) (RC3)
%
\bibitem[1989]{DDP89}
Disney M.J., Davies J.I., Phillipps S. 1989, MNRAS 239, 939 ({\bf DDP})
%
%\bibitem[1994]{DopRyd94}
%Dopita M.A., Ryder S.D. 1994, ApJ 430, 163
%
%\bibitem[1994]{Eva94}
%Evans R. 1994, MNRAS 266, 511
%
%\bibitem[1970]{Freeman}
%Freeman K.C. 1970, ApJ 160, 811
%
%\bibitem[1985]{Fro85}
%Frogel J.A. 1985, ApJ 298, 528
%
%\bibitem[1989]{Gil89}
%Gilmore G., King I., van der Kruit P.C. 1989, The Milky Way as a Galaxy
%(Geneva Observatory, Sauverny-Versoix, Switzerland)
%
%\bibitem[1994]{Gio94}
%Giovanelli R., Haynes M.P., Salzer J.J., Wegner G., Da Costa L.N.,
% Freudling W. 1994, AJ 107, 2036
%
%\bibitem[1994]{GouPhD}
%Goudfrooij P. 1994, Ph.D.\ Thesis, University of Amsterdam, The Netherlands
%
%\bibitem[1994]{Gou94}
%Goudfrooij P., Hansen L., J\o rgensen H.E., N\o rgaard-Nielsen H.U., de
%Jong T., van den Hoek L.B. 1994, A\&AS 104, 179
%
%\bibitem[1987]{Gre87}
%Green E.M., Demarque P., King C.R. 1987, The Revised Yale Isochrones and
%Luminosity Functions (Yale University Observatory, New Haven)
%
%\bibitem[1989]{GreAnd89}
%Grevesse N., Anders E. 1989, in: Cosmic Abundances of Matter, AIP
%Conference Proceedings 183, ed.\ C.J.~Waddington (AIP, New York)
%
%\bibitem[1992]{Han92}
%Han M. 1992, ApJS 81, 35
%
\bibitem[1941]{Henyey}
Henyey, L.G., Greenstein, J.L., 1941, ApJ, 93, 70
%
\bibitem[1994]{HuiPhD}
Huizinga J.E. 1994, Ph.D.~Thesis, University of Groningen, The Netherlands
%
\bibitem[1994]{Jan94}
Jansen R.A., Knapen J.H., Beckman J.E., Peletier R.F., Hes R. 1994, MNRAS 270, 373
%
%\bibitem[1996]{Jus96}
%Just A., Fuchs B., Wielen R. 1996, A\&A in press
%
%\bibitem[1989]{Ken89}
%Kennicutt Jr. R.C. 1989, ApJ 344, 685
%
%\bibitem[1984]{Kent84}
%Kent S.M. 1984, ApJS 56, 105
%
%\bibitem[1986]{Kent86}
%Kent S.M. 1986, AJ 91, 1301
%
%\bibitem[1987]{Kent87}
%Kent S.M. 1987, AJ 93, 816
%
\bibitem[1991]{Kna91}
Knapen J.H., Hes R., Beckman J.E., Peletier R.F. 1991, A\&A 241, 42
%
%%\bibitem[1995a]{Kna95a}
%%Knapen J.H., Beckman J.E., Heller C.H., Schlosman I., de Jong R.S. 1995a, ApJ
%%in press
%
%%\bibitem[1995b]{Kna95b}
%%Knapen J.H., Beckman J.E., Schlosman I., Peletier R.F., Heller C.H., 
%%de Jong R.S. 1995b, ApJ 443, L73
%
%%\bibitem[1986]{Kra86}
%%Kraan-Korteweg R.C. 1986, A\&AS 66, 255
%
\bibitem[1987]{KylBah87}
Kylafis N.D., Bahcall J.N. 1987, ApJ 317, 637
%
%\bibitem[1989]{ESO-LV}
%Lauberts A., Valentijn E.A. 1989, The Surface Photometry Catalogue of
%the ESO Sky Survey (ESO, Garching) (ESO-LV)
%
%\bibitem[1987]{LinPri87}
%Lin D.N.C., Pringle J.E. 1987, ApJ 320, L87
%
%\bibitem[1991]{MaeMey91}
%Maeder A., Meynet G. 1991, A\&AS 89, 451 
%
%\bibitem[1992]{Mat92}
%Mathewson D.S., Ford V.L., Buchhorn M. 1992, ApJS 81, 413
%
%\bibitem[1990]{Mat90}
%Mathis J.S. 1990, ARA\&A 28, 37
%
%\bibitem[1989]{Mat89}
%Matteucci F. 1989, in: Evolutionary Phenomena in Galaxies, eds.\ J.E.
%Beckman, B.E.J. Pagel (Cambridge University Press, Cambridge), p.~297
%
%\bibitem[1992]{Mat92}
%Matteucci F. 1992, in Morphological and Physical Classification of
%Galaxies, eds.\ G. Longo, M. Capaccioli, G. Busarello (Kluwer Academic
%Publishers, Dordrecht), p.~245
%
%\bibitem[1973]{Nilson} 
%Nilson P. 1973, Uppsala General Catalog of Galaxies   
%(Roy.~Soc.~Sci., Uppsala)  (UGC)
%
%\bibitem[1987]{OCon87}
%O'Connell R.W. 1987, in Stellar Populations, eds.\ C.~Norman, A.~Renzini,
%M.~Tosi (Cambridge University press, Cambridge), p.~167
%
%\bibitem[1989]{PelPhD}
%Peletier R.F. 1989, Ph.D.~Thesis, University of Groningen, The Netherlands
%
%\bibitem[1990a]{Pel90a}
%Peletier R.F., Davies R.L., Davis L.E., Illingworth G.D., Cawson M. 
%1990a, AJ 100, 1091
%
%%\bibitem[1990b]{Pel90b}
%%Peletier R.F., Valentijn E.A., Jameson R.F. 1990b, A\&A 233, 62
%
%\bibitem[1994]{Pel94}
%Peletier R.F., Valentijn E.A., Moorwood A.F.M., Freudling W. 1994,
%A\&AS 108, 621
%
%\bibitem[1995]{Pel95}
%Peletier R.F., Valentijn E.A., Moorwood A.F.M., Freudling W.,
%J.H. Knapen, J.E. Beckman 1995, A\&A Letters 300, L1
%
%\bibitem[1985]{Pic85}
%Pickles A.J. 1985, AJ 296, 340
%
%\bibitem[1986]{RenBuz86}
%Renzini A., Buzzoni A. 1986, in Spectral Evolution of Galaxies, eds.\ 
%C.~Chiosi, A.~Renzini (Reidel, Dordrecht)
%
\bibitem[1985]{RieLeb85}
Rieke G.H., Lebofsky M.J. 1985, ApJ 288, 618
%
%\bibitem[1993]{RixRie93}
%Rix H.-W., Rieke M.J. 1993, ApJ 418, 123
%
%\bibitem[1985]{Rub85}
%Rubin V.C., Burstein D., Ford Jr. W.K., Thonnard N. 1985, ApJ 289, 81
%
%\bibitem[1994]{RydDop94}
%Ryder S.D., Dopita M.A. 1994, ApJ 430, 142
%
%\bibitem[1955]{Sal55}
%Salpeter E.E. 1955, ApJ 121, 161
%
%%\bibitem[1993]{Sch93}
%%Schaerer D., Meynet G., Maeder A., Schaller G. 1993, A\&AS 98, 523 
%
%\bibitem[1976]{Sch76}
%Schweizer F. 1976, ApJS 31, 313
%
%\bibitem[1973]{Sea73}
%Searle L., Sargent W.L.W., Bagnuolo W.G. 1973, ApJ 179, 427
%
%\bibitem[1994]{SilEls94}
%Silva D.R., Elston R. 1994, ApJ 428, 511
%
%\bibitem[1990]{SomYos90}
%Sommer-Larsen J., Yoshii Y. 1990, MNRAS 243, 468
%
%\bibitem[1994]{Ter94}
%Terndrup D.M., Davies R.L., Frogel, J.A., DePoy D.L., Wells L.A. 1994,
%  ApJ 432, 518
%
%\bibitem[1980]{Tin80}
%Tinsley B.M. 1980, Fund.\ of Cos.\ Phys. 5, 287
%
%\bibitem[1977]{TulFis77}
%Tully R.B., Fisher J.R. 1977, A\&A 54, 661
%
%\bibitem[1990]{Val90}
%Valentijn E.A. 1990, Nat 346, 153
%
%\bibitem[1994]{Val94}
%Valentijn E.A. 1994, MNRAS 266, 614
%
%\bibitem[1985]{Alb85}
%van Albada T.S., Bahcall J.N., Begeman K.G., Sancisi R. 1985, ApJ 295, 305
%
%\bibitem[1986]{AlbSan86}
%van Albada T.S., Sancisi R. 1986, Phil.\ Trans.\ (Roy.\ Soc.\ London) A320, 447
%
%\bibitem[1985]{Van85}
%VandenBerg D.A. 1985, ApJS 58, 711
%
\bibitem[1988]{Kru88}
van der Kruit P.C. 1988, A\&A 192, 117
%
%\bibitem[1994]{AlvGer94}
%von Alvensleben U.A., Gerhard O.E. 1994, A\&A 285, 751
%
%\bibitem[1992]{VilEdm92}
%Villa-Costas M.B., Edmunds M.G. 1992, MNRAS 259, 121
%
\bibitem[1989]{Wai89}
Wainscoat R.J., Freeman K.C., Hyland A.R. 1989, ApJ 337, 163
%
%\bibitem[1986]{Wev86}
%Wevers B.M.H.R., van der Kruit P.C., Allen R.J. 1986, A\&AS 66, 505
%
\bibitem[1977]{Witt77}
Witt A.N., 1977, ApJS 35, 1
%
\bibitem[1992]{Witt92}
Witt A.N., Thronson Jr. H.A., Capuano Jr. J.M. 1992, ApJ 393, 611
%
%\bibitem[1994]{Wor94}
%Worthey G. 1994, ApJS 95, 107 ({\bf W94})
%
%\bibitem[1988]{WysGil88}
%Wyse R.F.G., Gilmore G 1988, AJ 95, 1404
%
%\bibitem[1989]{WysSil89}
%Wyse R.F.G., Silk J. 1989, ApJ 339, 700
%
%\bibitem[1994]{Zar94}
%Zaritsky D., Kennicutt Jr. R.C., Huchra J.P. 1994, ApJ 420, 87

\end{thebibliography}


\end{document}


