;+
; $Id: slabmodel2.pro,v 1.1 2009/04/22 18:32:58 antunes Exp $
;
; Project   : developer area for STEREO SECCHI
;                   
; Name      : slabmodel2
;               
; Purpose   : instantiates a slab model
;               
; Explanation: 
;               
; Use       : IDL> 
;    
; Inputs    : x,y,z,nel,modul,r,theta,phi,ii
;               
; Outputs   : returns the 3d model
;
; Keywords  : 
;               
; Calls from LASCO : none 
;
; Common    : none 
;               
; Restrictions: none
;               
; Side effects: 
;               
; Category    : models
;               
; Prev. Hist. : None.
;
; Written     : Paul Reiser, pre-2006
;               
; $Log: slabmodel2.pro,v $
; Revision 1.1  2009/04/22 18:32:58  antunes
; first half of relocation.
;
; Revision 1.1  2009/03/31 15:46:04  antunes
; More reorg
;
; Revision 1.1  2006/04/11 15:54:48  antunes
; Massive re-org of cvs 'dev' preparatory to moving into solarsoft.
;
; Revision 1.2  2006/01/09 17:14:42  antunes
; More commented added.
; New 'super3drv' driver routine for getting user specs for reconstruction
; process now checked in.  Just type 'super3drv' to test.
;
;-            
; -- density 14
; Sphericaly symmetric Ne
; use the Hayes model
; Ne(r)=Sum_k alpha_k * r^(-k) , k=1,2,3,4
; Orthoradial model : exp(-d^2/w0(r)) and exp(-abs(d/u0(r)))
; with w0(r) = Sum_k beta_k * r^(-k)), k=1,2,3,4
; see 11 May 2004
; the 5th parameter is a circular profile of 120 elements 
; used to modulate the density of the slab
; returns 1,2,3,4 if nel is zero, 0 if not

function slabmodel2,x,y,z,nel,modul,r,theta,phi,ii

;printf("\n check 4")

 DTOR=!pi/180.
 RADEG=1./DTOR

 x=float(x) 
 y=float(y)
 z=float(z)

 cvcoord,x,y,z,r,phi,theta

 if(r le 2.57)then begin
 nel=0.
 ret=1&goto,ret
 end

 if((phi lt 1.0472) or (phi gt 2.0944))then begin
 nel=0.
 ret=2&goto,ret
 end

; -- updated on June 15th using pminimizer02.pro
 c=[1.34e4,1.15e6,-6.022e6,5.577e7]; valid from 2.5 to 20 Rsun

 cw=[17.79*DTOR*DTOR,-85.93*DTOR*DTOR,138.32*DTOR*DTOR,391.45*DTOR*DTOR]; valid from 2.5 to 20 Rsun

; see 10 May 2004
 cu=[60.9441*DTOR , -542.687*DTOR , 1889.58*DTOR , -2289.55*DTOR]
 
 coef=[-1.0,-2,-3,-4]
 nel=0

; -- streamer half thickness
 w0=0.
 u0=0.
 
 for i=0,3 do begin
 ex=r^coef[i]
 nel=nel+c[i]*ex
 w0=w0+cw[i]*ex
 u0=u0+cu[i]*ex
 end

 dcross=w0/(2.*u0)
 
 if(abs(theta) lt dcross)then begin
 
 ee=(theta*theta)/w0
 if(ee gt 1E1)then begin
 nel=0.
 ret=3&goto,ret
 end
 
 nel=nel*exp(-ee) 
 
 end else begin

 ee=abs(theta/u0)-dcross/u0/2
 if(ee gt 1E1)then begin
 nel=0.
 ret=4&goto,ret
 end
 nel=nel*exp(-ee)
 end

 pos=fix(-119.+113.637*phi)

 if(pos lt 0)then pos=0
 if(pos gt 119)then pos=119

 nel=nel*modul[pos]

 ret=0

ret:

;print,x,y,z,r,theta,phi

 return, ret

end