pro rtcloud,plist,sbtot,imsize=imsize,$
fovpix=fovpix,$
obspos=obspos,obsang=obsang,$
nepos=nepos,neang=neang,$
c2image=c2image,c3image=c3image,$
quiet=quiet,$
cor1=cor1,cor2=cor2,hi1=hi1,hi2=hi2,c1fov=c1fov,$
hlonlat=hlonlat,$
instr=instr,rotmat=rotmat,$
obslonlat=obslonlat,$
projtype=projtype,rollang=rollang,$
scchead=scchead,pv2_1=pv2_1,pcin=pcin,flistout=flistout,$
fclip=fclip,listout=listout,listbehind=listbehind,$
nerotcntr=nerotcntr,nerotang=nerotang,nerotaxis=nerotaxis,$
netranslation=netranslation,unload=unload
;+
; $Id: rtcloud.pro,v 1.11 2013/01/29 17:50:58 thernis Exp $
;
; PURPOSE:
; Generate a view using a cloud of points
;
; CATEGORY:
; raytracing, simulation, 3d
;
; INPUTS:
; ! all the distances must be given in Rsun
; ! and the angles in rad
; plist : list of points : [3,N]
; imsize : [xs,ys] size of the output image
; 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
; rollang : allow to set the roll angle of the virtual instrument.
; Works only if a preset instrument is requested.
; -- 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
; nerotcntr : [x,y,z] center of rotation of the Ne model, in the Ne basis
; nerotang : [ax,ay,az] rotation of the Ne model around the nerotcntr, in the Ne basis
; nerotaxis : [axid1,axid2,axid3] axis id corresponding to the nerotang rotation angles. 1: X, 2: Y, 3: Z. Default is [3,2,1].
; netranslation : [tx,ty,tz] translation vector of the Ne model, in the Ne basis
; quiet : disable display of raytracing parameters
; 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
; projtype : projection type: (see Calabretta and Greisen,
; Representations of celestial coordinates in FITS, A&A
; 395, 1077-1122(2002))
; ARC : Zenithal equidistant (default)
; TAN : Gnomonic
; SIN : Slant orthographic
; AZP : Zenithal perspective
; If an instrument preset is requested then this keyword
; will overwrite the projection type of the selected
; instrument.
; pv2_1 : mu parameter for the AZP projection
; pcin : force the fits PCi_j matrix. Must be a 4 elements array
;
; instr : txt instrument preset, to select from the list above:
; scchead : secchi structure header: raytrace will use the
; positionning info of the header to generate the view
; flistout : set to 1 if you want a list of points instead of an image
; fclip : set to 1 if you want not to compute points masked by the sun
;
; -- Instrument FOV preset
; c1, c2, c3 : lasco C1, C2, C3
; cor1, cor2 : Secchi Cor1, Cor2
; hi1, hi2 : Secchi Hi1, Hi2
;
; OUTPUTS:
; sbtot : structure with image of the total brightness
; rotmat : final rotation matrix of the density cube
; listout : will contain the list of pixels if requested with /flistout
; listbehind : 1 if point masked by disk, else 0
;
;-
; -- use header information if passed by user
obslonlatflag=0L
if n_elements(scchead) ne 0 then begin
if scchead.instrume ne 'SECCHI' then flagsoho=scchead.telescop eq 'SOHO' else flagsoho=0b
wcs=fitshead2wcs(scchead,KEYWORD_NULL_VALUE=0)
if n_elements(instr) eq 0 then instr=scchead.detector
; if ~flagsoho && n_elements(secchiab) eq 0 then secchiab=strmid(scchead.obsrvtry,0,1,/reverse_offset)
if n_elements(obslonlat) eq 0 then begin
; ---- position of the virtual spacecraft
; crln_obs_rad = wcschangeunits(wcs.cunit[0],'rad',wcs.position.crln_obs)
; crlt_obs_rad = wcschangeunits(wcs.cunit[1],'rad',wcs.position.crlt_obs)
; obslonlat = float([crln_obs_rad, crlt_obs_rad, wcs.position.dsun_obs/(onersun()*1e3)])
obslonlat = float([(wcs.position.carr_earth + wcs.position.hgln_obs) * !dtor, wcs.position.crlt_obs * !dtor, wcs.position.dsun_obs / (onersun() * 1e3)])
; obslonlat = float([wcs.position.crln_obs*!dtor,wcs.position.crlt_obs*!dtor,wcs.position.dsun_obs/(onersun()*1e3)])
obslonlatflag=1L
endif
if n_elements(pv2_1) ne 0 then begin
pv2_1=float(scchead.pv2_1)
endif
if n_elements(rollang) eq 0 then begin
rollang=0.
endif
endif
; ---- define keyword default values
; and cast every inputs to avoid C crash with bad type entry
if n_elements(plist) eq 0 then plist=[0.,0.,0] else plist=float(plist)
if n_elements(imsize) eq 0 then imsize=[64L,64] else imsize=long(imsize)
if n_elements(fovpix) eq 0 then begin
fovpix=2./(64.)*!dtor
flagfovpix=0B
endif else begin
fovpix=float(fovpix)
flagfovpix=1B
endelse
if n_elements(obspos) eq 0 then begin
obspos=[0.,0,-214]
obsposflag=0b
endif else begin
obspos=float(obspos)
obsposflag=1b
endelse
if n_elements(obsang) eq 0 then begin
obsang=[0.,0,0]
obsangflag=0b
endif else begin
obsang=float(obsang)
obsangflag=1b
endelse
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(nerotcntr) eq 0 then nerotcntr=[0.,0,0] else nerotcntr=float(nerotcntr)
if n_elements(nerotang) eq 0 then nerotang=[0.,0,0] else nerotang=float(nerotang)
if n_elements(nerotaxis) eq 0 then nerotaxis=long([3,2,1]) else nerotaxis=long(nerotaxis)
if n_elements(netranslation) eq 0 then netranslation=[0.,0,0] else netranslation=float(netranslation)
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''',/info
; endelse
if n_elements(obslonlat) ne 0 and not obslonlatflag then begin
obslonlat=float(obslonlat)
obspos=obslonlat[2]*[sin(obslonlat[1]),$
sin(obslonlat[0])*cos(obslonlat[1]),$
-cos(obslonlat[0])*cos(obslonlat[1])]
obslonlatflag=1L
endif
if n_elements(obslonlat) eq 0 then begin
obslonlat=[-atan(obspos[1],obspos[2]),$
asin(obspos[1]/norm(obspos)),$
norm(obspos)]
obslonlatflag=0L
endif
if n_elements(rollang) eq 0 then rollang=0. else rollang=float(rollang)
if n_elements(flistout) eq 0 then flistout=0L else flistout=1L
if n_elements(fclip) eq 0 then fclip=0L else fclip=1L
; ---- detect instrument, if any
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
rtgetinstrwcsparam,instr,imsize,scchead,fovpix,crpix,obsangpreset,pc,projtypepreset=projtypepreset,pv2_1=pv2_1,rollang=rollang,crval=crval,pcin=pcin,flagfovpix=flagfovpix
; ---- compute incerse of pc matrix
pc=float(pc)
pcinv=float(invert(pc))
; ---- size of the list of points
if size(plist,/n_dim) eq 1 then nbp=1L else nbp=long((size(plist,/dim))[1])
listout=fltarr(2,nbp)
listbehind=lonarr(nbp)
; -- obsang and rollang is forced if passed by user
if not obsangflag and n_elements(instr) ne 0 then obsang=obsangpreset
if n_elements(pv2_1) eq 0 then pv2_1=0. else pv2_1=float(pv2_1)
; ---- set projection type
if n_elements(projtype) eq 0 then begin
if n_elements(projtypepreset) eq 0 then projtype='ARC' else projtype=projtypepreset
endif
projtype=strupcase(projtype)
case projtype of
'ARC' : projtypecode=1L
'TAN' : projtypecode=2L
'SIN' : projtypecode=3L
'AZP' : projtypecode=4L
else : message,'Bad projtype keyword !'
endcase
; ---- init the outputs
btot=fltarr(imsize[0],imsize[1])
rotmat=fltarr(3,3)
rtinitenv
starttime=systime(1)
s=call_external(getenv('RT_LIBFILE'),$
'rtcloud',$
imsize[0],imsize[1],$
fovpix,$
obspos,obsang,$
nepos,neang,$
btot,$
crpix,quiet,$
hlonlat,$
rotmat,obslonlat,obslonlatflag,projtypecode,pv2_1,pc,$
plist,nbp,pcinv,flistout,listout,fclip,listbehind,nerotcntr,$
nerotang,netranslation,nerotaxis,unload=unload)
if quiet eq 0 then begin
print,'Seconds ellapsed :'
print,systime(1)-starttime
endif
revision=getcvsrevision('$Revision: 1.11 $')
case 1 of
c2image : rtdetec='c2'
c3image : rtdetec='c3'
c1fov : rtdetec='c1'
cor1 : rtdetec='cor1'
cor2 : rtdetec='cor2'
hi1 : rtdetec='hi1'
hi2 : rtdetec='hi2'
else : rtdetec=''
endcase
sbtot={imsize:imsize,fovpix:fovpix,obspos:obspos,obsang:obsang,nepos:nepos,neang:neang,im:btot,hlonlat:hlonlat,revision:revision,rtdetec:rtdetec,rotmat:rotmat,obslonlat:obslonlat,projtype:projtype,pc:pc,pv2_1:pv2_1,wcs:wcs,nerotcntr:nerotcntr,nerotang:nerotang,netranslation:netranslation,nerotaxis:nerotaxis}
return
end
;
; CVSLOG:
; $Log: rtcloud.pro,v $
; Revision 1.11 2013/01/29 17:50:58 thernis
; Change computation of obslonlat variable.
;
; Revision 1.10 2011-08-09 21:49:30 thernis
; Remove unuseful keyword secchiab
;
; Revision 1.9 2011-08-03 14:40:49 thernis
; - cast pv2_1 to float to avoid problem with HI (AZP projection)
; - cast obslonlatflag to long
; - clean-up unneeded variables
;
; Revision 1.8 2010-08-26 13:23:59 mcnutt
; added KEYWORD_NULL_VALUE to call to fitshead2wcs
;
; Revision 1.7 2009/04/13 21:18:41 thernis
; Implement extra positioning parameters for the models.
;
; Revision 1.6 2009/03/06 22:14:42 thernis
; Implement neshift
;
; Revision 1.5 2008/09/23 14:16:20 thernis
; add unload kw
;
; Revision 1.4 2008/06/06 17:38:25 thernis
; Modif to make it work with SOHO data.
;
; Revision 1.3 2007/07/19 22:24:26 thernis
; Implement listbehind keyword
;
; Revision 1.2 2007/07/19 19:20:18 thernis
; Fix bug with obslonlat when it was forced by the user
;
; Revision 1.1 2007/07/10 21:12:19 thernis
; First commit.
;
;