;pro utm.pro
;universal transit modeller
;version 9/may/2003
;written by H.J. Deeg
;version history
;26/1/99 first working version
;27/1/99 losses are now calculated precisely, writing of lc to output file
;28/1/99 in gettranslum, the case 'if tscale ge lscale' is disabled. This  
        ;results in much better
	;precision for big occulters but ultimately it would be better to
	;re-enable this case and include loss-correction
;29/1/99 re-enabled above case and included correction calc for it.
;27/4/99 opt. reading of lightcurve (tflag); setting output-format (oflag)
;29/4/99 removed paramfile subroutines into separate files, changes to
         ;plotlc
;30/4/99 added input param 'tarr' and output param 'lumarr'
;15/11/99 changed case statement for tflag
;7/4/00  added input param 'epoch'
;11/4/00 added one digit precision to output lc
;10/7/00 fixed error that tflag=0 even if tarr is given in command-line
;10/7/00 still need to incorporate quad limbd from utmq 
;10/10/00 inserted obj_destroy at end to avoid leakage of obj vars
;5/2/01 inserted position angle 
;21/10/02 lchead param, allows prepending setupfile (header) to output lcs
;24/10/02 ooff param, adds offset to output lc
;9/5/2003 made unused mass parameter optional in setupfile

@setupfile ;compile this first

function str,num
;converts number to string and removes white space.
return,strcompress(string(num),/remove_all)
end


function readlc,filename,tflag
;reads file with time-points, alternatively (tflag = 2) read
;time - luminosity pairs (lightcurve)
datstr=strarr(50000)               ;can deal with 200 parameters
comstr=strarr(200)               ;lines with comments
ndat=0l & ncom=0 & tstr= ' '
openr,UN1,/get_lun,filename     ;do that here so file errors show up NOW
while(not eof(UN1)) do begin
    readf,UN1,tstr
    if (strmid(tstr,0,1) eq '#') or strlen(tstr) le 2 then begin
        comstr(ncom)=tstr
        ncom=ncom+1
   endif else begin
        datstr(ndat)=tstr
        ndat=ndat+1
    endelse
endwhile
if ndat ge 1 then datstr=datstr(0:ndat-1)
if ncom ge 1 then comstr=comstr(0:ncom-1)
;print,ndat,' parameter values read from file ',setfname
;print,ncom,' comment lines read'
free_lun,UN1

;define data arrays
tmarr=dblarr(ndat,tflag)
;read in the main data block
for k=0,ndat-1 do begin
    tok=str_sep(strtrim(strcompress(datstr(k)),2)," ")
; if n_elements(tok) ne tflag then print,"error in file ",filename," with line: ",tok

    tmarr(k,0)=double(tok(0))
    if tflag eq 2 then tmarr(k,1)=double(tok(1))
endfor
return,tmarr
end

pro plotxy,pos,maxdist,sunit
	maxpl=maxdist*1.05
   !x.range= [-maxpl,maxpl]
   !y.range= [-maxpl,maxpl]
;   !x.style=1  ;exact axis ranges
;   !y.style=1
   ;plot x-y plane
   wset,0 & plot,pos[*,0],pos[*,1],psym=1,xtitle='X ('+sunit+')',ytitle='Y ('+sunit+')'
   !y.range= [maxpl,-maxpl]  ;to put observer at bottom  
   ;plot x-z plane
   wset,4 & plot,pos[*,0],pos[*,2],psym=1,xtitle='X ('+sunit+')',ytitle='Z ('+sunit+')'
	oplot,[0.0],[!y.crange[0]],psym=4  ;symbol for observer's direction
end

pro plotlc,t,l,tunit,oflag
if max(t) ge 100000 then t=(t mod 100000) ;to avoid display problems with long numbers
!x.range= [min(t),max(t)]
wset,5

case oflag of
    0: begin
        lumstr='sys.lum'
        !y.range=[min(l),max(l)] 
        plot,t,l,xtitle='time ('+tunit+')',ytitle=lumstr,ystyle=2,xstyle=0,xtickformat='(g10.6)',xticks=3
    end
    1: begin
        lumstr='rel.loss'
        !y.range= [1,1e-6]
        plot,t,(l>1e-6),/ylog,xtitle='time ('+tunit+')',ytitle=lumstr,ystyle=1,xstyle=0,xtickformat='(g10.6)',xticks=3
    end
    2: begin
        lumstr='rel.mag'
        !y.range=[max(l),min(l)] 
        plot,t,l,xtitle='time ('+tunit+')',ytitle=lumstr,ystyle=2,xstyle=0,xtickformat='(g10.6)',xticks=3
    end
endcase

if oflag eq  1 then begin       ;rel loss in log scale
endif else begin                ;others in lin scale
endelse
end

;###DEFINE INIT METHODS
function body::init,tinit,num   ;inititalize body values
	;if keyword_set(size) then begin
	;self.radi=size
	;endif
	snum= str(num)
	self.num=num
	self.type=parasign(snum+"type")
        self.radi=float(parasign(snum+"radi"))
        self.mass=float(parasign(snum+"mass",1,/nowarn))
        self.period=double(parasign(snum+"period"))
        phaseyes=paratest(snum+"phase")
        epochyes=paratest(snum+"epoch")
        if epochyes then self.epoch=double(parasign(snum+"epoch"))
        if phaseyes and epochyes then begin
            print,"WARNING: epoch AND initial phase are given for object ",snum
            print," epoch will override phase"
            endif
        if epochyes then begin
            self.phase=(tinit-self.epoch)/self.period ;epoch overrides phase
            endif else begin  ;epoch not given - use phase and derive epoch
            self.phase=double(parasign(snum+"phase"))
            self.epoch=tinit-(self.period*self.phase)
            endelse
        self.inclin=double(parasign(snum+"inclin",90))
        self.rdist=double(parasign(snum+"rdist"))
        self.pa=double(parasign(snum+"pa",90))  ;position angle;dir: N towards E

	self -> geotcalc
	return, 1
end

;function planet::init,num
;	ival = self -> body::init(num)   ;call init 
;	snum= str(num)
;	self->print
;	return, 1
;end

function star::init,tinit,num
	snum= str(num)
	self.lum = float(parasign(snum+"lum"))
	self.limbd = float(parasign(snum+"limbd"))
;	self->print
	ival = self -> body::init(tinit,num)   ;call init 
	self ->geolcalc   ;calc luminosity distrib
	self.lformoc=float(self.lform)  ;copy init lum to occulted lum
	return, 1
end

function subbody::init,tinit,num
	snum= str(num)
	self.mainobj = fix(parasign(snum+"mainobj"))
	ival = self -> body::init(tinit,num)   ;call init body
	return, 1
end

function ring::init,tinit,num
	snum= str(num)
	self.transp = float(parasign(snum+"transp"))
	ival = self -> subbody::init(tinit,num)   ;call init subbody
	return, 1
end


;####DEFINE PRINT METHODS 
pro body::print   ;method to print body obj	
 	print,'# general values of object ',str(self.num)
;	print,'num',self.num
	print,'type ',self.type
	print,'radi',self.radi
	print,'mass',self.mass
	print,'period',self.period
	print,'phase',self.phase
	print,'inclination',self.inclin
        print,'pos angle',self.pa
	print,'rdist',self.rdist
	print,'positions',self.pos
end
pro star::print   ;method to print subclass moon
	self -> body::print
 	print,'# star specific'
	print,'lum',self.lum
	print,'limbd',self.limbd
end
pro subbody::print   ;method to print subclass moon
	self -> body::print
 	print,'# subbody specific'
	print,'mainobj',self.mainobj
end
pro ring::print   ;method to print subclass ring
	self -> subbody::print
 	print,'# ring specific'
;	print,'mainobj',self.mainobj
end


pro body::helpobj  ;prints help about this object and stops
	help,/obj,self
	stop
end

pro body::printpos
	print,format='("nr: ",i3," ",3f14.4," ",a)',self.num,self.pos,self.type
end


pro body::displformocinit  ;initializes display of lform array 
;        gsize=size(self.lformoc,/dimension) & gs=gsize[0] 
	window,self.num+10,xsize=(200),ysize=(200),xpos=0,ypos=(self.num*240),title=str(self.num)+' '+self.type+': luminosity'
end

pro body::disptforminit  ;initializes display of tform array 
;        gsize=size(self.lform,/dimension) & gs=gsize[0] 
	window,self.num+20,xsize=(200),ysize=(200),title=str(self.num)+' '+self.type+': transparency',xpos=220,ypos=(self.num*220)
end

pro star::displformoc  ;displays lformoc array 
	wset,self.num+10
	tvscl,0>(250<congrid(self.lformoc,200,200))
end

pro body::disptformoc  ;displays tform array 
	wset,self.num+20
	tvscl,0>(250<congrid(self.tform,200,200))
end


;###DEFINE FUNCTIONS RELATED TO POSITON
function subbody::getmainobjnum ;returns number of mainobj
	return,self.mainobj
end

function body::getpos   ;returns position of object
	return,self.pos
end

pro moon::subpos,mainpos 
        ;moon:adds mainpos to object's position-vector 
	self.pos=self.pos+mainpos
end
pro ring::subpos,mainpos 
        ;ring:takes mainpos as its own position
	self.pos=mainpos
end
	
pro body::incrempos,tinc,obj  ;sets/increments phase and postion 
    ;through timeinterval tinc. If tinc = 0.0, set posn at current phase
	if self.period ne 0.0 then self.phase=self.phase+tinc/self.period
	phi=self.phase*2*!pi
	tet=(90.0-self.inclin)*!pi/180.
        del=(90.0-self.pa)*!pi/180.

	xp=self.rdist*sin(phi)  	;posX'	
	yp=self.rdist*cos(phi)*sin(tet);posY'
	zp=self.rdist*cos(phi)*cos(tet);posZ'

        ;rotation matrix for position angle
	self.pos[0]=cos(del)*xp-sin(del)*yp  	;posX	
	self.pos[1]=sin(del)*xp+cos(del)*yp     ;posY
        self.pos[2]=zp          ;posZ


end

pro moon::incrempos,tinc,obj  ;position of moon is mainobj+moon pos
	self -> body::incrempos,tinc,obj
	mainobj=self.mainobj
	mainpos= obj(mainobj) -> getpos()   ;get pos of mainobj
	self.pos=self.pos+mainpos  	    ;add pos
end

pro ring::incrempos,tinc,obj     ;position of ring is mainobj pos!
	mainobj=self.mainobj
	mainpos= obj(mainobj) -> getpos()   ;get pos of mainobj
	self.pos=mainpos  	    ;update pos
end

;###DEFINE FUNCTION RELATED TO GEOM FORM
pro body::geotcalc
                                ;calculates bodies 2D tranparency distribution
gsize=size(self.tform,/dimension)  ;get size of form-array
gs3=gsize[0] * 3    		;array-size 3 times as big
gbig=bytarr(gs3,gs3)            ;make array 3 times as big
grad=gs3*0.45                    ;diameter will be ~90% of size of array
gctr=gs3/2                      ;x,y position-index of array center
grad2 = grad^2                  ;square of radius
for i=0,gs3-1 do begin          ;initialize gbig
    for j=0,gs3-1 do begin
        if ((i-gctr)^2+(j-gctr)^2) le grad2 then gbig(i,j) = 255 ;opaque
    endfor                      ;j
endfor                          ;i
                                ;the big array gbig is resized to gsize, 
				;and border pixels will be interpolated 
self.tform=rebin(gbig,gsize[0],gsize[1]) ;rebin to orig size
end

pro ring::geotcalc
                                ;calculates rings 2D tranparency distribution
gsize=size(self.tform,/dimension) ;get size of form-array
gs3=gsize[0] * 3                ;array-size 3 times as big
gbig=bytarr(gs3,gs3)            ;make array 3 times as big
grad=gs3*0.45                    ;outer diameter will be ~90% of size of array
girad=grad*self.rdist/self.radi
gctr=gs3/2 
grad2 = grad^2                  ;square of outer radius
girad2 = girad ^ 2              ;sq of inner radius
tet=(90.0-self.inclin)*!pi/180.
sintet=sin(tet)
transp=exp(alog(self.transp)/sintet) ;transparency in depend of inclin
transpc = fix(255*(1-transp))   ;here: 255: opaque, 0: transp
for i=0,gs3-1 do begin          ;initialize gbig
    for j=0,gs3-1 do begin
        r2= ((i-gctr)^2+((j-gctr)/sintet)^2)
        if r2  le grad2 and r2 ge girad2 then gbig(i,j) = transpc
    endfor                      ;j
endfor                          ;i
 ;the big array gbig is resized to gsize, and border pixels
 ;will be interpolated 
self.tform=rebin(gbig,gsize[0],gsize[1]) ;rebin to orig size
end

pro star::geolcalc
                                ;calculates bodies 2D luminosity distribution
gsize=size(self.tform,/dimension) ;get size of form-array
gs3=gsize[0] * 3		;array-size 3 times as big
;	gs4=(gs4- (gs4 mod 2))+1  ;increment gs4 to next odd number
gbig=bytarr(gs3,gs3)            ;make array 3 times as big
grad=gs3*0.45                    ;radius will be ~90% of size of array
gctr=gs3/2 
grad2 = grad^2                  ;square of radius
cl=self.limbd                   ;limbdarkening
for i=0,gs3-1 do begin          ;initialize gbig
    for j=0,gs3-1 do begin
        r2= ((i-gctr)^2+(j-gctr)^2) ;square of radius
        if r2 le grad2 then begin 
;	gbig(i,j) =255
            gbig(i,j) = byte(255*(1-cl*(1-sqrt(1-(r2/grad2))))) ;brightness distr
        endif
    endfor                      ;j
endfor                          ;i
                                ;the big array gbig is resized to gsize, and border pixels
                                ;will be interpolated 
self.lform=rebin(gbig,gsize[0],gsize[1]) ;rebin to orig size
          ;below conversion factor from pixels surface-brightness units(0-255)
          ;to luminosity unit (like Lsol)
self.lconv=self.lum/total(self.lform)
end

;### INFORMAL ROUTINES
function body::getrad   ;returns position of object
	return,self.radi
end

function body::getrdist   ;returns position of object
	return,self.rdist
end

function body::gettform   ;returns tform array of object
	return,self.tform
end

function star::getlum
	return,self.lum
end

;### RELATED TO OCCULTATION OF LIGHTPATH
function body::findinfront,i,zind,pos,rad
     ;return indices of objects in front of it
n_objects=n_elements(zind)
frind=-1   ;index of objects that are in front, -1: none
is=zind[i]                      ;sorted i index
for k=i+1,n_objects-1 do begin   ;from current object to closest in z
  ks=zind[k]                     ;sorted k index 	
  if abs(pos[is,0]-pos[ks,0]) lt rad[is] + rad[ks] then begin  ;test in x
    if abs(pos[is,1]-pos[ks,1]) lt rad[is] + rad[ks] then begin  ;test in y
      ;print,'object ',ks,' is in front of obj. ',is
      if frind[0] eq -1 then frind=[ks] else frind = [frind[*],ks]
    endif
  endif
endfor       
return,frind                  
end

function star::getlum  ;return intrinsic luminosity
       return, self.lum
end

pro star::resetlformoc   ;resets lformoc (occulted luminosity array)
       self.lformoc=float(self.lform)
end



function star::gettranslum,objfr ;calc luminosity with objects in array objfr in front

lxymax=dblarr(2) & lxymin=dblarr(2) ;min and max x-y extent of array of luminous body
lxymax=self.pos[0:1]+(self.radi*0.5/0.45)
lxymin=self.pos[0:1]-(self.radi*0.5/0.45)
txymax=dblarr(2) & txymin=dblarr(2) ;min and max x-y extent of transiting body
tindmin=intarr(2)               ;minimum index of subarray cutted out from tform 
tindmax=intarr(2)               ;maximum index  "   "
;tindctr=intarr(2)               ;center index  "
lindmin=intarr(2)               ;minimum index of subarray of lform corresponding to tform 
lindmax=intarr(2)               ;maximum index  "   "
lindctr=intarr(2)               ;center index  "

self.lformoc=float(self.lform)  ;reset lformoc array to intrinsic luminsos.
lsize=size(self.lform,/dimension) & lsize=lsize[0] ;get size of lform
lctr=lsize/2                    ;index of center pixel
lscale=self.radi/(0.45*lsize)   ;scale of lformoc array

n_object=n_elements(objfr)
for i=0,n_object-1 do begin
    lumorg = total(self.lformoc)   ;get luminosity of lformoc before transit  

    tpos=objfr[i] -> getpos()
    tform=objfr[i] -> gettform()
    trad=objfr[i] -> getrad()

    tsize=size(tform,/dimension) & tsize=tsize[0] ;get size of tform
    tctr = tsize /2             ;index of ctr pixel of tform
    tscale=trad/(0.45*tsize)
                                ;below the subarray indices are calculated
    tindmin=[0,0]> round((lxymin - tpos[0:1])/tscale + tctr)
    tindmax=[tsize-1,tsize-1] < round((lxymax-tpos[0:1])/tscale + tctr)
                                ;indices of tform within tindmin,tindmax are in front of self.lform
                                
                                ;corresponding indices of self.lform that correspond to tform
    lindctr= round((tpos[0:1]- self.pos[0:1])/lscale + lctr)
    lindradideal=trad*(0.5/0.45)/(lscale)   ;ideal size of subarray
    lindrad=(1>fix(lindradideal-0.5))  ;used size of subarray; - 0.5 accounts for ctr pix
    lindmin=[0,0]>(lindctr-lindrad)
    lindmax=[lsize-1,lsize-1] <(lindctr+lindrad)
;        lindmin=[0,0]> round((tpos[0:1]-trad - self.pos[0:1])/lscale + lctr)
;        lindmax=[lsize-1,lsize-1] < round((tpos[0:1]+trad-self.pos[0:1])/lscale + lctr)


;txymin=tpos[0:1]+((tindmin-tctr)*tscale);min and max xy positions of 
;txymax=tpos[0:1]+((tindmax-tctr)*tscale); tformoc

				;subarray of tform that is occulting lfrom
        tformoc=tform[tindmin[0]:tindmax[0],tindmin[1]:tindmax[1]]
                                ;subarray of lform that is behind occulting object 
        lformbh=self.lformoc[lindmin[0]:lindmax[0],lindmin[1]:lindmax[1]]
        lumbhorg=total(lformbh)  ;luminosity of lfrombh before transit

;				;below, depending on the scale of tform and lform,
;				;from each array (tform or lform) the overlapping 
;				;parts are scaled to
;				;the finer grid of the two.
;
    if tscale ge lscale then begin ;tform is coarser, and scaled to lform           
                                ;tformoc is the part of tform that is occulting lform
        tformocproj=congrid(tformoc,lindmax[0]-lindmin[0]+1,lindmax[1]-lindmin[1]+1)
                                ;apply transparency of occulter to lformbh
        lformbh = float(lformbh*float(255-tformocproj)/255.)
                                ;glue lformbh into complete lformoc array
	if total(lformbh gt 0.0) then begin ;do correction stuff
	lumbhafter=total(lformbh)
        ;the relative brightness loss that occured in the array lformbh
	rellosoccured=(lumbhorg-lumbhafter)/lumbhorg
        ;below is the scale in the projected array tformocproj 
	tscaleprojx=tscale*(tindmax[0]-tindmin[0]+1)/(lindmax[0]-lindmin[0]+1)
	tscaleprojy=tscale*(tindmax[1]-tindmin[1]+1)/(lindmax[1]-lindmin[1]+1)
	;tscaleprojx,y should be the same then lscale, but isn't due to integer array-size
	;For that reason, the loss is corrected by the ratio of these scales
	tscalecorrfac=(lscale^2)/(tscaleprojx*tscaleprojy)
        rellosdesired=rellosoccured/tscalecorrfac   ;that's the loss we want in lformbh
        ;below, lformbh is corrected to display rellosdesired   
        desiredlumbh=lumbhorg*(1-rellosdesired)
	lformbh=lformbh*(desiredlumbh/lumbhafter)
	endif   ;total(lformbh gt 0.0)
	;glue lfrombh back into lformoc
        self.lformoc[lindmin[0]:lindmax[0],lindmin[1]:lindmax[1]]= lformbh
;        print,'size tformoc',size(tformoc,/dimension) 
;        print,'size tformocproj',size(tformocproj,/dimension) 

    endif else begin 		;the lform array is scaled to tform ;this is the normal case
 	;scale lformbh to size of occulter's array
        lformbhproj=congrid(lformbh,tindmax[0]-tindmin[0]+1,tindmax[1]-tindmin[1]+1)
       	lumbhprojorg=total(lformbhproj)
	if lumbhprojorg gt 0.0 then begin  ;bypass of case 0.0 that can happen at begin or
                                           ;end of a transit 
        ;apply transparency of occulter to lformbhproj
        lformbhproj=float(lformbhproj*float(255-tformoc)/255.)
        lumbhprojafter=total(lformbhproj)
	;the relative brightness loss that occured in the projected array lformbhproj
	rellos=(lumbhprojorg-lumbhprojafter)/lumbhprojorg
	;below is the scale in the projected array lformbhproj 
	lscaleprojx=lscale*(lindmax[0]-lindmin[0]+1)/(tindmax[0]-tindmin[0]+1)
	lscaleprojy=lscale*(lindmax[1]-lindmin[1]+1)/(tindmax[1]-tindmin[1]+1)
	;lscaleprojx,y should be the same then tscale, but isn't due to integer array-size
	;For that reason, the loss is corrected by the ratio of these scales
	lscalecorrfac=(tscale^2)/(lscaleprojx*lscaleprojy)
        rellosdesired=rellos*lscalecorrfac   ;that's the loss we want in lformphproj and lformbh
        ;scale lformbhproj back to orig size of lformbh
        lformbh=congrid(lformbhproj,lindmax[0]-lindmin[0]+1,lindmax[1]-lindmin[1]+1,/minus_one)

	;the  lines below correct the value of lformbh to display the 
	;correct loss, which is relosscorr
        lumbhafter=total(lformbh)
	rellosoccured=(lumbhorg-lumbhafter)/lumbhorg
	desiredlumbh=lumbhorg*(1-rellosdesired)
	lformbh=lformbh*(desiredlumbh/lumbhafter)

	;glue lfrombh back into lformoc
        self.lformoc[lindmin[0]:lindmax[0],lindmin[1]:lindmax[1]]=lformbh

;the lines below are an other method to correct the loss, 
;that avoids the 'brighter square' problem, and
;that should be indentical to the used one.
;However, the errors in gsize=24 are notably bigger
;This may be due to calc. precision problem?
;	deltaoccured=(lumbhorg-lumbhafter)  ;the absolute loss (in pix-vals) that occured in lumbh
;	deltadesired=lumbhorg*rellosdesired ;the absolute loss that is desired
;	lumafteroccured=total(self.lformoc)   ;luminosity of lformoc after transit 
;	;calc corr factor to get lformoc to brightness it should have
;	corrfac=1D + double((deltaoccured - deltadesired)/double(lumafteroccured))
;	self.lformoc=self.lformoc*corrfac

	endif  ;if lumbhprojorg ge 0.0
     endelse	
endfor
;calc total luminosity
lumtot=total(self.lformoc)*self.lconv ; luminosity with occultation
return,lumtot
end

;###MAIN PROGRAM

pro utm,setfname,dflag,tarr,lumarr
common parrs,pvarr,pnarr

;read setupfile
readsetupfile,setfname

;get units
tunit=parasign("tunit")         ;get time-unit
;lunit=parasign("lunit")  ;get luminosity-unit
sunit=parasign("sunit")         ;get size-unit
;munit=parasign("munit")  ;get mass-unit


;initialize time-array and luminosity array
tflag = fix(parasign("tflag",0)) 
if n_elements(tarr) ne 0 then tflag =3 ;a time-array was given as a parameter
case 1 of
    tflag eq 0 :begin           ;make up times from tinit,tinc,tfin
;prototype routine for a clock
        tinc=double(parasign('tinc'))
        tinit=double(parasign('tinit'))
        tfin=double(parasign('tfin'))
        nlcpts=long((tfin-tinit)/tinc) +1
        tarr=dblarr(nlcpts)
        tcount = 0L & t=tinit
        for tcount=0L,nlcpts-1 do begin
            tarr(tcount)=t
            t=t+tinc
        endfor    
    end
    (tflag eq 1) or (tflag eq 2)  : begin ;tflag=1 or 2, read from file
        infilename=parasign('infilename')
        tmarr=readlc(infilename,tflag)
        tarr=tmarr[*,0]
        if tflag eq 2 then lum0arr=tmarr[*,1]
        nlcpts=n_elements(tarr)
    end
    tflag eq 3: nlcpts=n_elements(tarr) ;tarr is already defined as param
endcase
lumarr=fltarr(nlcpts)

;radfac=0.45  ;factor that describes ratio between size of object's
                                ;radius (in pixels) and corresponding array size 
;define structures for the object-types

gs=fix(parasign("gsize",12))    ;base-size of geometric arrays
gs=gs-(gs mod 2)                ;decrement to next even nr
;sizes of geometric arrays (tform,lform) to be used for various bodies 
;they will be odd numbers to have the center within a pixel:
gsplan=gs+1                     ;planet
gsring=gs*2+1                   ;ring
gsstar=gs*4+1                   ;star

a={body, num:0,type:'body',radi:0.,mass:0.,period:0D,phase:0D,epoch:0D,inclin:0D,pa:0D,rdist:0D,pos:dblarr(3)}
b={star,lum:0.0,limbd:0.0,lform:bytarr(gsstar,gsstar),lformoc:fltarr(gsstar,gsstar),lconv:0.0,tform:bytarr(gsstar,gsstar),inherits body}
c={planet,tform:bytarr(gsplan,gsplan),inherits body}
d={subbody,mainobj:0,inherits body}
e={moon,tform:bytarr(gsplan,gsplan),inherits subbody}
f={ring,transp:0.0,tform:bytarr(gsring,gsring),inherits subbody}

;create the objects and assign initial values to them

;first,find number of objects by looking for largest param Ntype
;that is defined in setupfile
n_object=0&while paratest(str(n_object)+"type") do n_object=n_object+1
if dflag ge 3 then print,str(n_object)+' objects will be defined'

;objects of type(i) are created and initialized
type=strarr(n_object)
obj=objarr(n_object)
for i=0,n_object-1 do begin
    type[i]=parasign(str(i)+"type")
    obj[i]=obj_new(type(i),tarr[0],i) ;call initialization method
    obj[i] -> incrempos,0.,obj  ;set initial position
    if obj_isa(obj[i],'star') then if dflag ge 3 then obj[i] -> displformocinit
endfor                          ;i=0,n_object


;CODE BELOW IS NOT WRITTEN OBJECT ORIENTED YET
if dflag ge 3 then window,0,xsize=400,ysize=320,title='x-y plane',xpos=400
if dflag ge 3 then window,4,xsize=400,ysize=320,title='x-z plane',xpos=800
if dflag ge 2 then window,5,xsize=800,ysize=320,title='light curve',xpos=400,ypos=320

;get some values of the objects that are needed in time-loop
zind=intarr(n_object)           ;sort-index in z-direction
pos=dblarr(n_object,3)          ;positions of all objects
rad=dblarr(n_object)            ;radii of all objects
for i=0,n_object-1 do rad[i] = obj[i] -> getrad()  
maxdist=0D
for i=0,n_object-1 do maxdist = (maxdist > obj[i] -> getrdist())  



;initialize display
for i=0,n_object-1  do if obj_isa(obj[i],'star') then if dflag ge 3 then obj[i]-> displformoc

;prepare output
oflag = fix(parasign("oflag",1))
case oflag of
    0: lumstr='sys.lum'
    1: lumstr='rel.loss'
    2: lumstr='rel.mag'
endcase
onoise=float(parasign("onoise",0,/nowarn)) ;opt noise to be added to output
ooff=float(parasign("ooff",0,/nowarn))  ;opt offset added to output
wflag = fix(parasign("wflag",0,/nowarn))
lchflag = strmid(strupcase(parasign("lchead","y",/nowarn)),0,1) eq "Y" ;check if setupfile (header) should be preprended to lc
if wflag eq 1 then begin
    lcoutname=parasign("lcfilename","lc.out")
    if lchflag then begin
        savesetupfile,lcoutname,/comm
        openw,unout,lcoutname,/get_lun,/append
        printf,unout,'#  '
        printf,unout,'# simulated lightcurve follows'
    endif else openw,unout,lcoutname,/get_lun
    printf,unout,'#   time ('+tunit+') '+ lumstr
endif



;get total system luminosity
syslum0 = 0.0
for i=0,n_object-1  do if obj_isa(obj[i],'star') then syslum0=syslum0+obj[i]->getlum()

lasttrflag=bytarr(n_object)     ;flags if there was a transit in last loop

;START MAIN LOOP
tcount=0L
for tcount=0L,nlcpts-1 do begin
    t=tarr[tcount]
    syslum=0.0
    for i=0,n_object-1 do pos[i,*] = obj[i] -> getpos()  
    zind=sort(pos[*,2])         ;contains sequence in z, farthest obj (-z) is first
    if dflag ge 4 then for i=0,n_object-1 do obj[i] -> printpos
    for i=0,n_object-1 do begin
        izs=zind[i]             ;current index within zind array
        if obj_isa(obj[izs],'star') then begin ;look for other objects in front of it
                                ; get array 'frind' with indexes of the objects that are in front
            frind= obj[izs] -> findinfront(i,zind,pos,rad)
                                ;frind contains either -1 (none in front) or the indexes
            if frind[0] eq -1 then begin  
                syslum=syslum+obj[izs] ->getlum() ;no transits 
                if lasttrflag[i] eq 1 then if dflag ge 3 then obj[izs] -> displformoc
                lasttrflag[i]=0
            endif else begin    ;transit case
                if dflag ge 4 then $
                  print,'in front of obj ',str(izs),' are obj(s) ',str(frind)
                syslum=syslum+ obj[izs] ->gettranslum(obj[frind]) 
                lasttrflag[i] = 1
                if dflag ge 3 then obj[izs] -> displformoc
                obj[izs] -> resetlformoc
            endelse
        endif
    endfor                      ;i
;    rloss=(syslum0-syslum)/syslum0

;   definitions of lumfin (output luminosities)    
    case oflag of
        0:lumfin=syslum         ;system luminosity
        1: lumfin=(syslum0-syslum)/syslum0 ;relative loss
        2: lumfin=(-2.5)*alog10(syslum/syslum0) ;magnitude changes
    endcase
    if onoise ne 0.0 then begin ;add noise to output
        lumfin=lumfin+randomu(seed,/normal)*onoise
    endif
    lumfin=lumfin+ooff  ;add offset
    if tflag eq 2 then lumfin=lumfin+lum0arr[tcount] ;add lumfin to input lc

    lumarr[tcount]=lumfin

    
    if dflag eq 1 then print,format='("t= ",f16.8," ",a,"= ",f9.6)',t,lumstr,lumfin
    if dflag ge 3 then plotxy,pos,maxdist,sunit
    if wflag then printf,unout,format='(f14.6," ",f9.6)',t,lumfin

    if dflag ge 2 then plotlc,tarr[0>(tcount-100):tcount],lumarr[0>(tcount-100):tcount],tunit,oflag 
    if dflag ge 4 then begin 
        tmp='' & read,tmp,prompt='hit return to continue'
    endif 
    if tcount lt nlcpts-2 then begin
        for i=0,n_object-1 do obj[i] -> incrempos,(tarr[tcount+1]-t),obj
    endif    
endfor                          ;t
if wflag then free_lun,unout
obj_destroy,obj                 ;kill object array
end