Calculating the dispersion of a solid core PCF

Hello,

In the last weeks I've been trying to calculate the dispersion (for the 
various modes) of a Blazephotonics NL-1.7-670 solid core PCF. I've read 
the documentation and older posts to this list and based on all that I 
made the ctl file attached that does a 2D calculation, with the k point 
in the z direction. The quantity k(omega) is obtained using the find-k 
function, and for each frequency an external program is called that 
creates the fiber profile in a .h5 file with the right permitivity so as 
to include the silica dispersion. The profile is then loaded using 
epsilon-input-file. To obtain the dispersion the quantity k(omega) is 
derivated using a 5 point finite difference formula. The problem is that 
the curve I get is the same as bulk silica. I made lots of changes to 
the fiber profile loaded, and couldn't get any of those to resemble the 
real fiber dispersion.

Is there anything wrong with this approach (ctl file attached)?

Addicionaly, if/when I calculate the dispersion for the various bands: 
are these bands TEM modes?

Thanks,
João Silva

_______________________________________________
mpb-discuss mailing list
[email protected]
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/mpb-discuss

(use-modules (goosh))
(define-param kz-min 0.01)
(define-param kz-max 40)
(set! num-bands 1)
(set! ensure-periodicity false)

; Length units: um
(set! geometry-lattice (make lattice (size 19.9355 19.2258 no-size)))

(set! resolution 32)

(define k_values 0.0)

;------------------------------------------------------------------------
;Calculates the square of silica's refraction index using a Sellmeier
;formula
;------------------------------------------------------------------------
(define (Silica_eps l)	; l - Wavelength in um
	(+ 1 
	    (/ 0.6961663 (- 1.0 (sqr (/ 0.0684043 l)))) 	
	    (/ 0.4079426 (- 1.0 (sqr (/ 0.1162414 l)))) 	
	    (/ 0.8974794 (- 1.0 (sqr (/ 9.896161 l)))) 	
	)
)

(define (loop_omega x x-max dx)
   (if (<= x x-max)
      (begin
	(run_prog "./pnm_to_h5" (number->string (Silica_eps (/ 6.28318530718 x))))
	(set! epsilon-input-file "nl17_mask.h5")


	(init-params NO-PARITY true)
	 (set! k_values (find-k NO-PARITY x 1 num-bands (vector3 0 0 1) 1.0E-4 1.8 kz-min kz-max ))
         (print "dispk\t" x "\t" (first k_values) "\n")
         (loop_omega (+ x dx) x-max dx))))

(loop_omega 4.1 12.5 0.1) ; execute loop from a to b in steps of dx