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.
;
;