NRT home work #4

Implementing a formal RT solver

You can choose one of the three algorithms: Feautrier, DELO or Hermitean. The necessary data can be found in the following binary file.

nrt_hw4.bin is a FORTRAN binary file. FORTRAN binary file consists of one or several records. Each record has the following structure: <length><data><length> where length is a 4-byte integer representation of the record length in bytes.

nrt_hw4.bin contains the following data:

Record1: Ndepth - (integer*4) number of depth points;

Nwave - (integer*4) number of wavelength points;
STEP - (real*8) stepsize in depth in centimeters;

Record2: WAVE - (real*8) array of Nwave elements, wavelengths;

Record3: OPAC - (real*4) array of Nwave elements, opacities;
SRC - (real*4) array of Nwave elements, source function;
INIT_INT - (real*4) array of Nwave elements, boundary cond;

Record4 (repeats Ndepth-1 times):
OPAC - (real*4) array of Nwave elements, opacities;
SRC - (real*4) array of Nwave elements, source function

Reading nrt_hw4.bin in FORTRAN*77 may look like this:

INTEGER*4 Nd,Nw,Id,Iw

REAL*8 STEP,WAVE(600)

REAL*4 OPAC(600),SRC(600),INIT_INT(600)

. . .

OPEN(UNIT=1,FILE='nrt_hw4.bin',FORM='UNFORMATTED',STATUS='OLD')

READ(1) Nd,Nw,STEP

READ(1) (WAVE(Iw),Iw=1,Nw)

READ(1) (OPAC(Iw),Iw=1,Nw),(SRC(Iw),Iw=1,Nw),(INIT_INT(Iw),Iw=1,Nw)

DO Id=2,Nd

. . .

READ(1) (OPAC(Iw),Iw=1,Nw),(SRC(Iw),Iw=1,Nw)

. . .

END DO

CLOSE(1)

. . .

FORTRAN*90 version only differs in the declaration section:

INTEGER(KIND=4) :: Nd,Nw,Id,Iw

REAL(KIND=8) :: STEP,WAVE(600)

REAL(KIND=4) :: OPAC(600),SRC(600),INIT_INT(600)

. . .

OPEN(UNIT=1,FILE='nrt_hw4.bin',FORM='UNFORMATTED',STATUS='OLD')

READ(1) Nd,Nw,STEP

READ(1) (WAVE(Iw),Iw=1,Nw)

READ(1) (OPAC(Iw),Iw=1,Nw),(SRC(Iw),Iw=1,Nw),(INIT_INT(Iw),Iw=1,Nw)

DO Id=2,Nd

. . .

READ(1) (OPAC(Iw),Iw=1,Nw),(SRC(Iw),Iw=1,Nw)

. . .

END DO

CLOSE(1)

. . .

In IDL one can do more things:

. . .

openr,1,'nrt_hw4.bin',/F77

Nd=0L

Nw=0L

STEP=0.D0

READU,1,Nd,Nw,STEP

WAVE=dblarr(Nw)

READU,1,WAVE

OPAC=fltarr(Nw)

SRC=fltarr(Nw)

INIT_INT=fltarr(Nw)

READU,1,OPAC,SRC,INIT_INT

FOR Id=1,Nd-1 do begin

. . .

readu,1,OPAC,SRC

. . .

ENDFOR

close,1

. . .

Here is something to compare with

I include my DELO calculations so you can compare. My results are in FORTRAN binary file which contains 3 records:

Record #1: Nwl (integer*4) number of wavelengths
Record #2: WL (real*8) wavelengths, array of Nwl elements
Record #3: SP (real*4) synthetic spectrum, array of Nwl elements

The corresponding observed solar spectrum in a similar format can be found here

Deadline for presenting/handing in your results is November 26th 2012!