1. Introduction

In a photonic-bandgap fiber (PBF), the cladding is a photonic crystal made of a two-dimensional periodic refractive index structure that creates bandgaps where propagation at certain optical frequencies is forbidden. This structure can be either a series of concentric rings of alternating high and low indices, as in so-called Bragg fibers, [1] or a two-dimensional periodic lattice of air holes arranged on a geometric pattern (e.g., triangular), as in air-core PBFs. [2–6] A defect is introduced in the periodic structure in the form of a larger air hole (the core), which introduces guided modes within the previously forbidden bandgap. PBFs exhibit a number of unique optical properties, including nonlinearities about two orders of magnitude lower than in a conventional fiber, a greater control over the mode properties, and interesting dispersion characteristics, such as modal phase velocities greater than the speed of light in air. For these reasons, in optical communications and sensor applications alike PBFs offer exciting possibilities previously not available from conventional fibers.

Since the propagation of light in air-core fibers is based on a very different principle than in a conventional fiber, it is important to develop a sound understanding of their still poorly understood properties, to improve these properties, and to tailor them to specific applications. Since their fairly recent invention,[2–6] modeling has been an important tool for analyzing the modal properties of air-core PBFs, particularly their phase and group velocity dispersion,[7–9] transmission spectra,[4,6] unique surface modes,[10–14] and propagation loss.[15–17] However, extensive and systematic parametric studies of these fibers are lacking because simulations of these complex structures are conducted with analytical methods that are either time-consuming or not applicable to all PBF geometries. For example, simulating all the modes of a typically single-mode air-core PBF with the commonly used MIT Photonic Bandgap (MPB) code, [18, 19] with a cell size of 10Λ × 10Λ and a spatial resolution of Λ/16 (Λ is the period of the photonic-crystal cladding) currently takes approximately 10 hours on a supercomputer using 16 parallel 2-GHz processors. Faster computations have been reported by White et al using a multipole decomposition method, with which they modeled all modes in a hexagonal PBF with four rings of holes at one wavelength in only about one hour on a 733-MHz personal computer. [20] This figure would be reduced to a few minutes on a faster PC (e.g., ~15 minutes on a 3.2-GHz processor). Poladian et al demonstrated a polar-coordinate decomposition method that modeled the fundamental mode (effective index and fields) in 6 minutes.[21] However, this high speed was accomplished by limiting the number of air hole layers to one, and it still required a supercomputer. Finite-difference and finite-element methods have also been developed, [22, 23] and they are quite fast and versatile, as described further on. Also, much finer resolutions are needed to accurately model the extremely small physical features of air-core fibers, in particular their very thin membranes.

In the light of these limitations, it is of critical importance to develop computer codes that give the ability to model PBFs of any geometry (and other waveguides with complex structures and fine features) quickly and accurately. In this paper, we describe a new mode solver that meets this goal by applying the finite-difference method of Refs. [22] and [23] and the shift-invert technique, described in Ref. 24, to solve a magnetic field eigenvalue equation in a sparse matrix form to find the modes (core, surface, ring, and bulk modes) of a PBF of arbitrary index profile. This simulator, which we call the Stanford Photonic Bandgap Fiber (SPBF) code, also incorporates a novel technique to reduce the bandwidth of the sparse matrix approach, which further reduces computation time by at least a factor of 10. Compared to a plane-wave expansion, this approach reduces computation time by about three orders of magnitude and it cuts down the amount of required data storage by about two orders of magnitude, so that full simulations of all the mode properties of a given fiber can now be carried out very quickly on a personal computer. This method is also simple to program. The mode solver program, implemented on a PC, is only 200 lines of MATLAB code.

2. State of the art

Photonic-bandgap fibers are typically modeled using numerical methods originally developed for the study of photonic crystals. Most methods are based on solving the photonic-crystal master equation:

where r is a vector (x, y) that represents the coordinates of a particular point in a plane perpendicular to the fiber axis, ε(r) = n ²(r) is the dielectric permeability of the fiber cross-section at this point, where n(r) is the two-dimensional refractive index profile of the fiber, H(r) is the mode’s magnetic field vector, ω is the optical frequency, and c the speed of light in vacuum. This formulation is most useful for 3D photonic-bandgap structures, and it is solved by assuming a constant value for the propagation constant k_z and computing the eigenfrequencies ω for this value. Six basic methods have been used to solve for modes in the context of air-core PBFs. The first one is the MIT Photonic Bandgap software (MPB). [18, 19] The solutions of Eq. (1) are obtained in the spatial Fourier domain, so the first step is to decompose the mode fields and 1/ε(r) in spatial harmonics. Equation (1) is then written in a matrix form and solved by finding the eigenvalues of this matrix. This is done by setting k_z to some value k₀ and solving for the frequencies at which a mode with this wave-number occurs. The second method is the beam-propagation method (BPM). [25] A pseudo-random field is launched into the fiber, which excites different modes with different effective indices. The propagation axis is transformed into an imaginary axis, so that the effective indices become purely imaginary and the fiber modes become purely attenuating. Instead of accumulating phase at different rates, the modes now attenuate at different rates, which give us the ability to separate them and extract them iteratively by propagating them through a length of fiber. The mode with the lowest effective index, which attenuates at the smallest rate, is extracted first, and the other modes are computed in order of increasing effective index. The third method is the finite difference time domain method (FDTD), which solves the time-dependent Maxwell’s equations:

Solving these equations requires considering the time variable as well. [26] The fourth method involves using a two-dimensional mode eigenvalue equation, where the eigenvalue is the mode effective index. The mode is then expanded in plane-waves or in wavelet functions [24, 27], the refractive index profile is decomposed in Fourier components, and the resulting matrix equation is solved for eigenvalues and eigenvectors. The multipole method [20] yields extremely accurate and quick calculations of the modes of PBFs with circular or elliptical hole shapes. However, it is limited to these shapes; it does not work, for example, for PBFs with scalloped cores or near-hexagonal cladding holes. The last two methods are the finite-difference [22, 23] and finite-element methods, where the waveguide index profile is digitized and the gradients involved in the mode equation approximated by finite-difference equations. These methods have been successfully applied to holey fibers and other waveguides, and they offer the speediest calculation.

When applied to the calculation of the modes of a PBF, all these methods suffer from a number of disadvantages. First, with the exception of the multipole, finite-difference, and finite-element methods, computation times can be very long. A second disadvantage is that techniques that rely on a mode expansion require taking the Fourier (or wavelet) transform of the fiber refractive index profile. Since this profile contains very fine features (thin membranes) and abrupt discontinuities (at the air-glass interfaces), both of which are key to the mode behavior, a large number of harmonics should ideally be retained to faithfully describe these features. Because this would require much more memory and computer processing power than is available, even in supercomputers, the discontinuities in the index profile are smoothened, which artificially widens the membranes. This standard practice results in significant differences between the index profile that is being modeled and the actual profile, and it introduces systematic errors in the predictions. A third limitation is that some existing mode solvers must calculate the modes of the highest order bands first and work their way down, which propagates computation errors. This requirement also increases computation time, since it forces the computer to calculate more modes than the ones of interest. The multipole method is the most powerful of the existing methods, but unfortunately it only works speedily for a few specific hole shapes, circular and elliptical mostly, thus making it less versatile.

3. The SPBF method

The objective of this work is to model the fields and propagation constants of the modes (core, cladding, surface modes, etc.) of an optical waveguide with an arbitrary index profile n(x, y) that is translation invariant along the fiber longitudinal axis (z). The transverse magnetic field h _T(r) of the waveguide modes satisfies the eigenvalue equation: [28]

where $k_0=\frac{2\pi }{{\lambda }_0}$(λ₀ is the wavelength), h _T = h_x(x, y) ux + h_y(x, y) uy is the transverse magnetic field, h_x and h_y are the components of the transverse magnetic field projected onto the unit vectors u_x and u_y , ${\nabla }_T=\frac{\partial }{\partial x}{\mathbf{u}}_x+\frac{\partial }{\partial y}{\mathbf{u}}_y$ and n_eff is the mode effective index. With this notation, the total magnetic field vector H defined earlier is simply h _T + h_z u_z , where u_z is the unit vector along the fiber z-axis. Equation (3) is satisfied by all modes (TE, TM, and hybrid modes). The longitudinal component of the magnetic field and the components of the electric field can all be calculated from the transverse magnetic field using:

Equation (3) is vectorial and has two equations that couple the two components h_x (x, y) and h_y (x, y) of the transverse magnetic field, and thus they must be solved simultaneously. It is important to note that the form of Eq. (3) is such that it can be solved it by fixing the wavelength and calculating all the fiber modes at this wavelength.

By using a finite-difference technique based on Yee’s cell to account for the field and refractive index discontinuities, as well as index-averaging over each discretization pixel that straddles an air-silica boundary, as done in [29], using the entire set of Maxwell’s equations, it is possible to discretize Eq. (3) with a high accuracy. The use of index averaging also allows improving the accuracy of the computation, at the cost of smoothing the refractive index structure over a small percentage of the total number of pixels. [29]

After discretization, Eq. (3) can thus be written as:

where Π is a linear operator acting on a vector field. To solve this equation, we discretize the eigenvalue problem by sampling the two fields h_x (x, y) and h_y (x, y) in a plane cross-section perpendicular to the fiber axis on a grid of m_x points in the x direction and m_y points in the y direction. This operation yields a new magnetic field vector h(x, y) of dimension 1 × 2m_xm_y , which contains (arranged using some arbitrary numbering system) S = m_x m_y samples of h_x (x, y) and S samples of h_y (x, y). This sampling is done over an area that respects the crystal periodicity, i.e., over a minimum cell that is a fundamental component of the cladding’s photonic-crystal structure. Similarly, the linear operator Π is sampled using finite-difference equations (which entails, in particular, sampling the refractive index profile n (x, y) of the fiber), so that it is now represented by a matrix M of dimension 2S × 2S. Equation (5) is solved subject to standard boundary conditions. Unlike with a k-domain method, which relies on the Bloch theorem and thus assumes periodic boundary conditions, with this new method one of three types of boundary conditions can be used: (1) fields vanishing outside the simulation area, (2) fields periodically repeated over the entire space, with the simulation domain as a unit super cell, or (3) surrounding the simulation domain by a perfectly matched layer (for propagation loss calculation). Both the boundary conditions and the simulation domain boundary are taken into account when sampling the operatorΠ, so that this information is also contained in the matrix M itself. Equation (4) is thus replaced by a single matrix equation:

Finding the roots of Eq. (6), i.e., the eigenvalues (n_eff ) and eigenmodes (field distributions) of the fiber modes, is a typical matrix eigenvalue problem. However, the operator matrix M is typically extremely large. For example, consider the case of a photonic-bandgap fiber with a triangular pattern of circular air holes of radius ρ = 0.47Λ, and a circular core of radius R = 0.8Λ, 4 rows of air holes in the cladding, and a square cladding boundary. The cell is then sixteen periods on the side, i.e., its dimension is 10Λ × 10Λ. To resolve the fiber’s thin membranes and achieve sufficient accuracy, we might want to use a grid step size (or resolution) of Λ/50. The number of sampling points is then m_x = 10 × 50 in x and m_y = 10 × 50 in y. Thus S = m_x m_y = 250,000, and M will contain (2S) ² = 250 billion elements. Finding the eigenvalues of matrices this large is out of the range of most computers, as even storing them is problematic. However, the replacement of the differential operators by their finite-difference equivalents makes M a sparse matrix. This is illustrated in Fig. 1, which shows the general form of the matrix representing the mode equation for such an example. The horizontal and vertical axes represent the coordinate of the elements in the matrix, which each go up to a maximum value of 2S = 516,000. The solid dots, which merge into solid curves because their density is too high to be individually resolved on this diagram, represent nonzero elements. The blanks between them are all zero elements. In this particular example, the total number of non-zero elements n_z is 2,796,826.

 

Fig. 1. Example of a matrix M representing the mode equation for a particular PBF using a rectangular boundary.

Download Full Size | PDF

The density D of this kind of matrix, defined as the ratio of non-zero elements to the total number of elements is $\frac{5}{2S}\le D\le \frac{7}{2S},$ the exact value depending on the index profile. In the above example, the matrix density is only of the order of D ≈ 1.05 × 10^-5, or equivalently only about 0.001% of the matrix elements are non-zero. The implication is that in spite of the large size of the matrix, the number of elements that must be stored during computation is comparatively low, down from ~266 billion to ~2.8 million for this typical example, which is well within the range of what can be handled with standard personal computers.

Another critical parameter in the matrix M is its bandwidth B_M , defined as the maximum distance between a non-zero coefficient and the first diagonal.

A small bandwidth is imperative for the fast computation of eigenvalues in sparse matrices. Since the bandwidth of the original matrix, M is proportional to S and quite large, as illustrated in Fig. 1, it is desirable to lower it. By applying a proper permutation [30] on the elements of both h and M, it is possible to rearrange M so that all the non-zero elements are clustered near the first diagonal and the bandwidth is minimized. This process yields a new sparse matrix N and a new permutated magnetic field vector g, which also satisfy Eq. (6), i.e., N g = $k_0^2$ $n_{\mathit{\text{eff}}}^2$. The bandwidth of N is now proportional to S ^1/2 and greatly reduced from the original value. This is illustrated in Fig. 2, which shows the matrix N obtained by applying this process to the original matrix M of Fig. 1. Most of the non-zero elements are now close to the diagonal. The bandwidth has been reduced from 460,320 (Fig. 1) to only 1,922 (Fig. 2), i.e., by a factor of ~240, which speeds up the calculation of the eigenvalues of the matrix by about an order of magnitude.

 

Fig. 2. Matrix N obtained by rearranging matrix M of Fig. 1 to minimize the matrix bandwidth.

Download Full Size | PDF

In principle, we could find the modes of the fiber by diagonalizing N to find all of its eigenvalues. However, this approach would require calculating and storing a very large number of eigenvalues and eigenvectors, which would take a great deal of time and memory. Furthermore, this method is generally wasteful because it calculates all of the eigenvalues. Instead, a more expedient way to solve Eq. (6) is to calculate only the few eigenvalues of interest, as done for example in Ref [24]. This restriction is not unreasonable, since in general one is primarily interested in finding only a few modes within the bandgap. The problem can thus be re-stated as calculating the m eigenvalues of N that are closest to a given effective index value n ₀, where m is a small number. For example, to model the bandgap edges and effective indices of the modes within the bandgap of a particular fiber at a particular wavelength, we might be interested in calculating the m ₀ = 10 modes that are closest to n ₀ = 1, since the core modes of an air-core fiber have effective indices close to 1. The next issue is that calculating the eigenvalues of a large matrix that are closest to a given value is still a difficult task. To overcome this problem, as in Ref. [24] we make use of the property that the m ₀ eigenvalues of a matrix N that are closest to some value n ₀ are the same as the largest (in amplitude) eigenvalues of the matrix:[24]

where I is the identity matrix. This property is useful here because the eigenvalues of maximum amplitude of a given matrix can be calculated very quickly by using what is known as the Courant-Fischer theorem.[31] This theorem restricts its eigenvalue search to small blocks of dimension m ₀ within the matrix, which is much faster than identifying all the eigenvalues. This approach is applied by first calculating a new matrix N - $k_0^2$ $n_0^2$ I , then inverting it using a standard LU decomposition (the matrix is expressed as the product of a lower and an upper triangular matrices) to obtain matrix A, which is fast because N - $k_0^2$ $n_0^2$ I is also a sparse matrix with a minimum bandwidth, and finally making use of the Courant-Fischer algorithm to calculate only the m ₀ largest eigenvalues of A, which are the effective indices of interest. We have therefore transformed our problem from looking for a specific set of eigenvalues around a given value, which is a difficult problem, to looking for the largest eigenvalues of a transformed matrix, which is considerably faster. [24]

 

Fig. 3. Summary of the SPBF code algorithm (see text for details).

Download Full Size | PDF

The steps in the mode-solving algorithm outlined above are summarized in Fig. [3]. The first step is to generate the sparse matrix M representing the action of Maxwell’s equations on the transverse magnetic field. This requires using the sampled refractive index profile distribution of the fiber as well as the wavelength and the boundary conditions. In the next step, M is rearranged to yield a new sparse matrix N with a minimum bandwidth to speed up the extraction of the eigenvalues. The third step is to specify the effective index n ₀ around which we want to calculate the PBF modes. We then translate the spectrum of N by calculating N-$k_0^2$ $n_0^2$ I. The fourth step is to invert this matrix through LU decomposition. The final step is to compute the largest eigenvalues and eigenvectors of this matrix using the Courant-Fischer theorem.

This mode-solving algorithm can be implemented in a straightforward manner using any one of a number of commercial mathematical software programs. For example, the mode solver we describe in the following treatment was written in MATLAB. Most of the required mathematical operations are readily available from MATLAB, in particular a library named UMFPACK that computes the largest eigenvalues of a matrix using the Courant-Fischer theorem mentioned above. Thanks to these built-in features, and to the relative simplicity of the basic mode-solving algorithm, only minimal programming is required, and the overall program contains less than 200 lines of code.

4. Features and performance of the SPBF code

The advantages of this new mode solver are multiple. First, it allows considerably speedier and easier simulations of PBFs than other simulators, and it is does not require a supercomputer. When run on a personal computer with a 3.2-GHz Pentium IV processor and 4 gigabytes of RAM, the SPBF code finds the propagation constants and fields of 20 modes in a given fiber in around 4 minutes per wavelength for a rectangular grid 550 points on the side. Treating the same problem using the MPB code [18, 19] takes around 1 hour using 16 2-GHz parallel processors, and it cannot be done in a reasonable time on a personal computer. Our code uses the wavelength as the variable parameter, a more intuitive approach traditionally used to model conventional fibers and optical waveguides. It also calculates only the required modes, and it can be used with a variety of boundary conditions. It can model refractive index distributions with negative or complex values, which makes it possible to estimate mode propagation losses using perfectly matching boundary conditions. Another useful feature is that if one needs to zoom in on a particular detail of a mode’s field distribution, for example the tail in one of the air holes surrounding the core, one can use the results of the calculation using a coarser grid as new boundary conditions to extend the calculation to a finer grid, as is possible with other finite-difference and finite-element methods. This process should give access to very fine features of PBF modes.

5. Simulation of air-core and holey fibers

 

Fig. 4. Sampling of the refractive index profile of an air-core PBF with an air-filling ratio ρ = 0.47Λ and a core radius R = 0 8Λ used to simulate this fiber using periodic boundary conditions.

Download Full Size | PDF

We illustrate the performance of the SPBF code by applying it to model the propagation constants and intensity profiles of selected modes of an air-core PBF with circular air holes in a triangular lattice in silica. The fiber cross-section used in the simulation is shown in Fig. 4. The air-hole radius is ρ = 0.47Λ and the radius of the circular air core is R = 0.8Λ. The gray areas represent air and the black areas silica. The modeling was done using periodic boundary conditions: the simulation domain (super cell) represented in Fig. 4 was repeated periodically to tile the entire space and the modes of this periodic waveguide were calculated. The separation distance between two fiber cores in this periodic tiling was selected large enough (10Λ) to ensure negligible interaction between neighboring fiber cores and thus accurate predictions. This means that the simulation domain is a rectangle 10Λ wide (a little under 5 rings), a size commonly used in simulations [14] because it offers a good compromise between accuracy and computing time and memory requirements. The sampling of the fiber refractive index distribution, and thus the digitization of the fiber features (circular core and cladding air holes) was done with a spatial resolution of Λ/50. This finite resolution accounts for the small irregularities visible at the edges of the holes in Fig. 4.

 

Fig. 5. Mode intensity profile of the fundamental core mode of the PBF of Fig. 4 at λ = 0.60Λ. The contour lines are magnified in the right-side figure for easier viewing.

Download Full Size | PDF

 

Fig. 6. Cut at x = 0 of the fundamental core mode intensity of Fig. 5.

Download Full Size | PDF

Figure 5 shows the contour map of the intensity profile of one of the fundamental modes, for a first PBF with a core radius R = 0.8Λ. As is well known, the core modes are mainly concentrated in air, and this mode, being degenerate with an orthogonally polarized mode, does not exhibit a perfect six-fold symmetry. The intensity profile (cut along x = 0) is shown in Fig. 6. This profile exhibits lobes at the air-silica core boundaries (y = ±1.5Λ, see Fig. 6), as pointed out by other authors.

The dispersion curve of a PBF identical to that of Fig. 4, except with a core radius R = 1.0Λ, is shown in Fig. 7. The bandgap extends approximately from λ = 0.57Λ to λ = 0.645Λ. This bandgap supports two fundamental core modes, and no surface modes, a result consistent with the predictions of [13, 14] obtained using MPB.

 

Fig. 7. Band profile and dispersion of the two fundamental modes of the PBF of Fig. 4, except for a core radius R = 1.00Λ.

Download Full Size | PDF

 

Fig. 8. Convergence of the fundamental mode effective index vs. grid size for the PBF in Fig. 4; the reference value used is the result obtained by SPBF with 600 grid points per side.

Download Full Size | PDF

 

Fig. 9. Convergence of the fundamental mode effective index vs. grid size for a silica fiber in air (red circles) and for the holey microstructured optical fiber (MOF) of Ref. 22 (blue dots); the analytical solution is used as a reference for the silica fiber, and the multipole result as a reference for the MOF.

Download Full Size | PDF

To demonstrate the convergence of the SPBF code, we first show in Fig. 8 how the calculated value of the effective index converges as the number of grid points is increased for the fundamental core mode of the PBF in Fig. 4: the quantity plotted is |n_eff - n_eff |, where the reference value n_ref is the effective index calculated by SPBF with 600 grid points per side. The code reaches convergence to the fourth decimal place for a step size of Λ/42 (~420 points per side), and fifth decimal place for a step size of Λ/50 (~500 points per side).

To assess the accuracy of the SPBF code, we used it to calculate the effective index of the fundamental mode of: (1) a silica core of 3 μm radius, surrounded by air, at a wavelength of 1.5 μm, and compared to the reference value obtained by means of the analytical solution; (2) the holey microstructured optical fiber (MOF) of Ref. 22, illustrated in Fig. 3 of that reference, and compared to the reference value calculated in Ref. [22] with a high accuracy using the multipole method. The result of these comparisons is plotted in Fig. 9, where the effective index values calculated with the SPBF code are shown as red circles for the fiber, and blue dots for the MOF. The convergences rates as a function of the grid size are very similar in both cases, reflecting the effectiveness of the averaging technique. For grid sizes in the range of 100 to 200 points per axis, the SPBF code and the correct answers agree to below the fifth decimal place in both cases, and reach an agreement to within less than 10^-7 for grid sizes exceeding 450 points. This result confirms both the accuracy and the convergence of the SPBF code. It should be noted that the SPBF code achieves a higher accuracy (by two orders of magnitude) for each of these two fibers than for the PBF of Fig. 4 because they have a higher mode confinement, which enabled us to reduce the simulation window size, and thus the step size, which in turn improves the code’s accuracy.

 

Fig. 10. Band profile and dispersion of the two fundamental modes and surface modes of the same PBF as shown in Fig. 4 except for a core radius R = 1.15Λ.

Download Full Size | PDF

Figure 10 shows the calculated dispersion diagram of a second PBF, identical to the previous one, shown in Fig. 4, except that the core radius is increased to R = 1.15Λ. This fiber exhibits surface modes, which is consistent with previously published results. [13, 14] As shown in the intensity contour map of Fig. 11, these surfaces modes are localized at the interface between the air core and first layer of cladding air holes. Their fields are maximum in silica and evanescent in air. The surface mode shown in Fig. 11 also exhibits the six-fold symmetry expected from the fiber symmetry. These results are consistent with previously reported PBF mode behavior. [10–12]

 

Fig. 11 Intensity profile of an exemplary surface mode of the fiber modeled in Fig. 8 (ρ = 0.47Λ and R = 1.15Λ) at a wavelength λ = 0 607Λ The contour lines are magnified in the right-side figure for easier viewing.

Download Full Size | PDF

For the purpose of simulation of PBFs, the best way to minimize the total computation time using SPBF is to solve for a large number of modes of the structure (around 40) around a target effective index n ₀ = 1 or 0.99, until a core mode is found. This first step is necessary since we do not know the effective index of core modes prior to simulation. Once the core mode effective index has been determined at a given wavelength, this effective index, or an effective index calculated from it by linear extrapolation, can be used as the target effective index at the next wavelength. The number of modes solved at this wavelength can thus be lowered, to typically less than 20, since in general only the core modes, surface modes, and a few bulk modes at the edges of the bandgap need to be determined. This procedure reduces the computation time by tracking a particular mode of interest, and it centers the calculated modes on the core modes.

6. Conclusions

We have developed a new numerical method that allows us to simulate the modes of air-core photonic-bandgap fibers of arbitrary index profiles extremely quickly on a PC, while reducing the amount of data storage required by about two orders of magnitude. The effective indices and field profiles of the modes of a typical PBF at a particular wavelength can now be calculated in only 15 seconds per mode on a personal computer with 4 gigabytes of RAM. These substantial improvements were accomplished by (1) solving a vectorial transverse-magnetic-field equation in a matrix form, which can be done quickly because this matrix is sparse; (2) re-arranging the matrix elements to reduce the matrix bandwidth, an original technique which speeds up eigenvalue finding by another factor of at least one order of magnitude, and (3) solving exclusively for the modes of interest. This method is also simple to program (about 200 lines of MATLAB code), it can be used in conjunction with one of several kinds of boundary conditions, and it is applicable to quickly calculate the modes of any optical fiber or waveguide, regardless of the complexity of its index profile. The accuracy achieved for the mode effective index is better than the fifth decimal for a large structure like an air-core fiber, and better than the seventh decimal for smaller structures such as microstructured fibers. We hope to make the Stanford Photonic-Bandgap Fiber code available soon to scientists and researchers either on Stanford University’s website or through Stanford's Office of Technology Licensing.