; ---- wrapper that calls the C raytracing program pro rtwlseg,sbtot,sbpol,snetot,imsize=imsize,$ fovpix=fovpix,segment=segment,angle=angle,nbpix=nbpix,$ obspos=obspos,obsang=obsang,$ nepos=nepos,neang=neang,$ losnbp=losnbp,losrange=losrange,$ modelid=modelid,modparam=modparam,$ save=save,file=file,fakelasco=fakelasco,$ c2image=c2image,c3image=c3image,quiet=quiet,$ cor1=cor1,cor2=cor2,hi1=hi1,hi2=hi2,c1fov=c1fov,$ unload=unload,hlonlat=hlonlat,xdr=xdr,secchiab=secchiab,$ instr=instr,limbdark=limbdark,obslonlat=obslonlat ;+ ; $Id: rtwlseg.pro,v 1.3 2006/11/09 19:43:30 thernis Exp $ ; ; PURPOSE: ; Perform raytracing along a straight line segment profile ; ; CATEGORY: ; raytracing, simulation, 3d ; ; INPUTS: ; ! all the distances must be given in Rsun ; ! and the angles in rad ; segment : [[x1,y1],[x2,y2]] defines the segment extremities, ; in image pix. ; angle : if defined the segment goes from the image center and has ; nbpix lenght. Angle is overread by segment. ; nbpix : size of the output profile in pix ; imsize : [xs,ys] size of the image (not calculated here, just for ; the scale) ; fovpix : fov angle of one pixel in rad ; -- observer position and attitude ; obspos : [x,y,z] position of the observer in the Sun basis ; obslonlat : [lon,lat,height] position of the observer in Carrington ; coordinate. If set, then obspos is ignored. The optical ; axis always points toward the Sun center. Use obsang to ; change telescope orientation. ; obsang : [ax,ay,az] orientation of the observer, ; z is the optical axis ; -- Ne position and attitude ; nepos : [x,y,z] position of the Ne reference in the Sun basis ; neang : [ax,ay,az] orientation of the Ne ; -- LOS params ; losnbp : number of step for the integration along the LOS ; losrange : [lstart,lend] range for the integration along the LOS ; in Rsun. The origin of the LOS is the orthogonal ; projection of the Sun cntr on that LOS. ; modelid : model id number ; modparam : parameters of the model ; save : put path and filename in that variable (without extention) ; if you want to save the results in a .fits binary table. ; file : density cube filename for density 4 and 5 ; fakelasco : put fake lasco header information in the fits header ; quiet : disable display of raytracing parameters ; unload : set to force unloading the C++ program (actualy .so lib): ; useful when recompiling the C program ; hlonlat : [Hlon,Hlat,Hrot] heliographic lon and lat of the center of ; the disk, rotation angle corresponding to the projection of the ; north pole, counterclockwise. ; xdr : save into xdr format instead of fits table. 'save' keyword ; must be set for xdr to take effect. ; secchiab : 'A' or 'B', to select Ahead or Behind spacecraft, for ; secchi only ; instr : txt instrument preset, to select from the list above: ; -- Instrument FOV preset ; c1, c2, c3 : lasco C1, C2, C3 ; cor1, cor2 : Secchi Cor1, Cor2 ; hi1, hi2 : Secchi Hi1, Hi2 ; limbdark : limb darkening coeff: default 0.58 ; ; OUTPUTS: ; sbtot : structure for the total brightness profile ; sbpol : structure for the polarized brightness profile ; snetot : structure for the LOS integrated electron density profile ; ;- ; ---- define keyword default values ; and cast every inputs to avoid C crash with bad type entry if n_elements(imsize) eq 0 then imsize=[64L,64] else imsize=long(imsize) if n_elements(segment) ne 0 then segment=float(segment) if n_elements(angle) ne 0 then angle=float(angle) if n_elements(nbpix) eq 0 then nbpix=255L else nbpix=long(nbpix) if n_elements(fovpix) eq 0 then fovpix=2./(64.)*!dtor else fovpix=float(fovpix) if n_elements(obspos) eq 0 then obspos=[0.,0,-214] else obspos=float(obspos) if n_elements(obsang) eq 0 then obsang=[0.,0,0] else obsang=float(obsang) if n_elements(nepos) eq 0 then nepos=[0.,0,0] else nepos=float(nepos) if n_elements(neang) eq 0 then neang=[0.,0,0] else neang=float(neang) if n_elements(losnbp) eq 0 then losnbp=64L else losnbp=long(losnbp) if n_elements(losrange) eq 0 then losrange=[-3.2,3.2] else losrange=float(losrange) if n_elements(modelid) eq 0 then modelid=1L else modelid=long(modelid) if n_elements(modparam) eq 0 then modparam=0. else modparam=float(modparam) if n_elements(quiet) eq 0 then quiet=0L else quiet=1L if n_elements(hlonlat) eq 0 then hlonlat=[0.,0,0] else hlonlat=float(hlonlat) if n_elements(secchiab) eq 0 then secchiab='A' else begin secchiab=strupcase(secchiab) if secchiab ne 'A' and secchiab ne 'B' then message,'secchiab keyword must be either ''A'' or ''B''' endelse if n_elements(limbdark) eq 0 then limbdark=0.58 else limbdark=float(limbdark) if n_elements(obslonlat) ne 0 then begin obslonlat=float(obslonlat) obslonlatflag=1L endif else begin obslonlat=[0.,0,0] obslonlatflag=0L endelse c2image=keyword_set(c2image) if c2image then instr='c2' c3image=keyword_set(c3image) if c3image then instr='c3' c1fov=keyword_set(c1fov) if c1fov then instr='c1' cor1=keyword_set(cor1) if cor1 then instr='cor1' cor2=keyword_set(cor2) if cor2 then instr='cor2' hi1=keyword_set(hi1) if hi1 then instr='hi1' hi2=keyword_set(hi2) if hi2 then instr='hi2' xdr=keyword_set(xdr) ; ---- instrument presets if requested rtinstrpresets,instr,imsize,fovpix,obsang,imszratio,secchiab=secchiab ; ---- define the default segment if n_elements(segment) eq 0 and n_elements(angle) ne 0 then begin center=rtcalccntr(fovpix,obsang,imsize) point=center+float(nbpix)*[cos(angle),sin(angle)] segment=[[center],[point]] endif else segment=[[0.,0.],[10,0]] btot=fltarr(nbpix) bpol=fltarr(nbpix) netot=fltarr(nbpix) crpix=fltarr(2) case modelid of 4 : begin print,'Model 4:' print,'Loading YM density...' print,'Raytracing will be done with trilinear intepolation' load_dens,dens,rco,phico,thetaco,file=file modparam=[rco,phico,thetaco,reform(dens,n_elements(dens))] dens=0b end 5 : begin print,'Model 5:' print,'Loading YM density...' print,'Raytracing will be done with nearest neighbor' load_dens,dens,rco,phico,thetaco,file=file modparam=[rco,phico,thetaco,reform(dens,n_elements(dens))] dens=0b end else : endcase rtinitenv starttime=systime(1) s=call_external(getenv('RT_LIBFILE'),$ 'rtwlseg',$ imsize[0],imsize[1],$ fovpix,$ obspos,obsang,$ nepos,neang,$ losnbp,losrange,modelid,$ btot,bpol,netot,modparam,$ crpix,nbpix,segment,quiet,hlonlat,limbdark,$ obslonlat,obslonlatflag,unload=unload) if quiet eq 0 then begin print,'Seconds ellapsed :' print,systime(1)-starttime endif ; -- sun radius in pix solar_r=abs(atan(1./obspos[2]))/fovpix mkhdr,h0,btot sxaddpar,h0,'COMMENT','Generated with raytrace software - NRL' ; -- X simu -> Y image, Y simu -> X image sxaddpar,h0,'CRPIX1',crpix[1]*solar_r+(float(imsize[0])-1)/2.,'X Sun center (pix)' sxaddpar,h0,'CRPIX2',crpix[0]*solar_r+(float(imsize[1])-1)/2.,'Y Sun center (pix)' sxaddpar,h0,'SOLAR_R',solar_r,'pix' sxaddpar,h0,'CDELT1',fovpix*!radeg*3600.,'arcsec' sxaddpar,h0,'CDELT2',fovpix*!radeg*3600.,'arcsec' sxaddpar,h0,'SOLAR_B0',0. sxaddpar,h0,'FOVPIX',fovpix,'rad' case 1 of c2image : rtdetec='c2image' c3image : rtdetec='c3image' c1fov : rtdetec='c1fov' cor1 : rtdetec='cor1' cor2 : rtdetec='cor2' hi1 : rtdetec='hi1' hi2 : rtdetec='hi2' else : rtdetec=0 endcase sxaddpar,h0,'RTDETEC',rtdetec,'raytrace detector option' ; -- put fake date ;sxaddpar,h0,'TIME_OBS','00:00:00.000' ;sxaddpar,h0,'DATE_OBS','2004/01/22' ; -- put fake LASCO header if keyword_set(fakelasco) then begin sxaddpar,h0,'FILEORIG','040124_023058.img' sxaddpar,h0,'DATE','2004/02/19 04:38:31.577' sxaddpar,h0,'DATE-OBS','2004/01/24' sxaddpar,h0,'TIME-OBS','02:30:05.451' sxaddpar,h0,'P1COL',20 sxaddpar,h0,'P1ROW',1 sxaddpar,h0,'P2COL',1043 sxaddpar,h0,'P2ROW',1024 sxaddpar,h0,'VERSION',2 sxaddpar,h0,'EXPTIME',1. sxaddpar,h0,'EXP0',23.000000 sxaddpar,h0,'EXPCMD',23.000000 sxaddpar,h0,'EXP1',1.8842800 sxaddpar,h0,'EXP2',2.0659200 sxaddpar,h0,'EXP3',2.5966800 sxaddpar,h0,'TELESCOP','SOHO' sxaddpar,h0,'INSTRUME','LASCO' sxaddpar,h0,'DETECTOR','C2' if keyword_set(c3image) then sxaddpar,h0,'DETECTOR','C3' sxaddpar,h0,'READPORT','C' sxaddpar,h0,'SUMROW',0 sxaddpar,h0,'SUMCOL',0 sxaddpar,h0,'LEBXSUM',1 sxaddpar,h0,'LEBYSUM',1 if imsize[0] eq 512 then begin sxaddpar,h0,'LEBXSUM',2 sxaddpar,h0,'LEBYSUM',2 endif if imsize[0] eq 256 then begin sxaddpar,h0,'LEBXSUM',4 sxaddpar,h0,'LEBYSUM',4 endif if imsize[0] eq 128 then begin sxaddpar,h0,'LEBXSUM',8 sxaddpar,h0,'LEBYSUM',8 endif sxaddpar,h0,'SHUTTR',0 sxaddpar,h0,'LAMP',0 sxaddpar,h0,'FILTER','Orange' if keyword_set(c3image) then sxaddpar,h0,'FILTER','Clear' sxaddpar,h0,'POLAR','Clear' sxaddpar,h0,'LP_NUM','Normal' sxaddpar,h0,'OS_NUM',3389 sxaddpar,h0,'IMGCTR',433 sxaddpar,h0,'IMGSEQ',0 sxaddpar,h0,'COMPRSSN','X%' sxaddpar,h0,'HCOMP_SF',64 sxaddpar,h0,'FP_WL_UP',0.0000000 sxaddpar,h0,'FP_WL_CM',0.0000000 sxaddpar,h0,'WAVELENG',0.0000000 sxaddpar,h0,'FP_ORDER',0 sxaddpar,h0,'M1_PZ1',0 sxaddpar,h0,'M1_PZ2',0 sxaddpar,h0,'M1_PZ3',0 sxaddpar,h0,'MID_DATE',53028 sxaddpar,h0,'MID_TIME',9018.2500 sxaddpar,h0,'PLATESCL',11.900000 if keyword_set(c3image) then sxaddpar,h0,'PLATESCL',56. sxaddpar,h0,'OFFSET',600 sxaddpar,h0,'IMAGE_CT',433 sxaddpar,h0,'SEQ_NUM',0 sxaddpar,h0,'OBT_TIME',1.4536026e+09 sxaddpar,h0,'R1COL',20 sxaddpar,h0,'R1ROW',1 sxaddpar,h0,'R2COL',1043 sxaddpar,h0,'R2ROW',1024 sxaddpar,h0,'BUNIT',' 0' sxaddpar,h0,'EFFPORT','C' sxaddpar,h0,'RECTIFY','TRUE' sxaddpar,h0,'CRVAL1',0.0000000 sxaddpar,h0,'CRVAL2',0.0000000 sxaddpar,h0,'CROTA',0.0000000 sxaddpar,h0,'XCEN',0.0000000 sxaddpar,h0,'YCEN',0.0000000 sxaddpar,h0,'CROTA1',180.00000 sxaddpar,h0,'CROTA2',180.00000 sxaddpar,h0,'CTYPE1','SOLAR-X' sxaddpar,h0,'CTYPE2','SOLAR-Y' sxaddpar,h0,'CUNIT1','ARCSEC' sxaddpar,h0,'CUNIT2','ARCSEC' sxaddpar,h0,'SECTOR',' 0' sxaddpar,h0,'RSUN',0.0000000 sxaddpar,h0,'NMISSING',0 sxaddpar,h0,'MISSLIST',' 0' endif sxaddpar,h0,'SEGMENTX1',segment[0,0],'pix' sxaddpar,h0,'SEGMENTY1',segment[1,0],'pix' sxaddpar,h0,'SEGMENTX2',segment[0,1],'pix' sxaddpar,h0,'SEGMENTY2',segment[1,1],'pix' sxaddpar,h0,'OBSPOSX',obspos[0],'Rsun' sxaddpar,h0,'OBSPOSY',obspos[1],'Rsun' sxaddpar,h0,'OBSPOSZ',obspos[2],'Rsun' sxaddpar,h0,'OBSANGX',obsang[0],'rad' sxaddpar,h0,'OBSANGY',obsang[1],'rad' sxaddpar,h0,'OBSANGZ',obsang[2],'rad' sxaddpar,h0,'NEPOSX',nepos[0],'Rsun' sxaddpar,h0,'NEPOSY',nepos[1],'Rsun' sxaddpar,h0,'NEPOSZ',nepos[2],'Rsun' sxaddpar,h0,'NEANGX',neang[0],'rad' sxaddpar,h0,'NEANGY',neang[1],'rad' sxaddpar,h0,'NEANGZ',neang[2],'rad' sxaddpar,h0,'LOSNBP',losnbp,'points' sxaddpar,h0,'LOSRNG1',losrange[0],'Rsun' sxaddpar,h0,'LOSRNG2',losrange[1],'Rsun' sxaddpar,h0,'MODELID',modelid if n_elements(file) ne 0 then $ sxaddpar,h0,'DENSFILE',file else $ sxaddpar,h0,'DENSFILE','' sxaddpar,h0,'HLONLAT0',hlonlat[0],'rad' sxaddpar,h0,'HLONLAT1',hlonlat[1],'rad' sxaddpar,h0,'HLONLAT2',hlonlat[2],'rad' if cor1 or cor2 or hi1 or hi2 then sxaddpar,h0,'SECCHIAB',secchiab,'Spacecraft A or B' hbtot=h0 sxaddpar,hbtot,'UNITS','Bsun' sxaddpar,hbtot,'COMMENT','Total brightness' hbpol=h0 sxaddpar,hbpol,'UNITS','Bsun' sxaddpar,hbpol,'COMMENT','Polarized brightness' hnetot=h0 sxaddpar,hnetot,'UNITS','e- . cm-3' sxaddpar,hnetot,'COMMENT','Integrated electron density' revision=getcvsrevision('$Revision: 1.3 $') sbtot={imsize:imsize,nbpix:nbpix,segment:segment,fovpix:fovpix,obspos:obspos,obsang:obsang,nepos:nepos,neang:neang,losnbp:losnbp,losrange:losrange,modelid:modelid,crpixprg:crpix,im:btot,fitshead:hbtot,hlonlat:hlonlat,rtdetec:rtdetec,densfile:(n_elements(file) ne 0 ? file : ''),revision:revision,secchiab:secchiab} sbpol={imsize:imsize,nbpix:nbpix,segment:segment,fovpix:fovpix,obspos:obspos,obsang:obsang,nepos:nepos,neang:neang,losnbp:losnbp,losrange:losrange,modelid:modelid,crpixprg:crpix,im:bpol,fitshead:hbpol,hlonlat:hlonlat,rtdetec:rtdetec,densfile:(n_elements(file) ne 0 ? file : ''),revision:revision,secchiab:secchiab} snetot={imsize:imsize,nbpix:nbpix,segment:segment,fovpix:fovpix,obspos:obspos,obsang:obsang,nepos:nepos,neang:neang,losnbp:losnbp,losrange:losrange,modelid:modelid,crpixprg:crpix,im:netot,fitshead:hnetot,hlonlat:hlonlat,rtdetec:rtdetec,densfile:(n_elements(file) ne 0 ? file : ''),revision:revision,secchiab:secchiab} if n_elements(save) ne 0 then begin print,'Saving the results' if xdr then begin savestruct,save+'btot.xdr',sbtot savestruct,save+'bpol.xdr',sbpol savestruct,save+'netot.xdr',snetot endif else begin mwrfits,sbtot,save+'btot.fits',/create mwrfits,sbpol,save+'bpol.fits',/create mwrfits,snetot,save+'netot.fits',/create endelse endif return end ; ; CVSLOG: ; $Log: rtwlseg.pro,v $ ; Revision 1.3 2006/11/09 19:43:30 thernis ; Implement obslonlat keyword to allow positioning the observer using Carrington spherical coordinates. ; ; Revision 1.2 2006/10/30 21:52:40 thernis ; Change the environment variable that points to the shared object library. ; ;