pro surf__define
; Store surface parameters
;
; CMR August 2003
;
; units
;          irho - radial variable for derivatives (1=rhon, 2=an, 3=psin)
;           rho - effective radius defined as sqrt(toroidal flux/pi B0) (m)
;        rhomax - maximum effective radius                              (m)
;        drhodr - d/dr (rho)
;          rhon - normalised effective radius
;             a - local minor radius                                    (m)
;          amax - maximum minor radius                                  (m)
;          dadr - d/dr (a) 
;            an - normalised minor radius
;           psi - poloidal flux psi-psi_core                            (Wb)
;        psimax - range of poloidal flux from edge to core              (Wb)
;        dpsidr - d/dr (psi)                                            (Wb)
;          psin - normalised poloidal flux
;         kappa - elongation
;      dkappadr - d/dr (kappa)
;         delta - triangularity
;      ddeltadr - d/dr (delta)
;             q - safety factor
;          dqdr - d/dr (q)
;             p - pressure                                              (Pa)
;          dpdr - d/dr (p)                                              (Pa)
;        d2pdr2 - (d/dr)^2 (p)                                          (Pa)
;         omega - toroidal rotation frequency                        (rad/s)
;      domegadr - d/dr omega                                         (rad/s)
;          rhoi - ion Larmor radius                                      (m)
;          rhoe - electron Larmor radius                                 (m)
;            ni - ion density                                         (m^-3)
;          dene - electron density                                    (m^-3)
;            ti - ion temperature                                       (eV)
;            te - electron temperature                                  (eV)
;          etai - L_ni/L_Te
;          etae - L_ne/L_Te
;            b0 - magnetic field                                         (T)
   dummy={surf,irho:0L,rho:-9.9d9,rhomax:-9.9d9,drhodr:-9.9d9,rhon:-9.9d9,a:-9.9e9,amax:-9.9e9,dadr:-9.9e9,an:-9.9e9,psi:-9.9e9,psimax:-9.9e9,dpsidr:-9.9e9,psin:-9.9e9,kappa:-9.9e9,dkappadr:-9.9e9,delta:-9.9e9,ddeltadr:-9.9e9,q:-9.9e9,dqdr:-9.9e9,p:-9.9e9,dpdr:-9.9e9,d2pdr2:-9.9e9,omega:-9.9e9,domegadr:-9.9e9,rhoi:-9.9e9,rhoe:-9.9e9,ni:-9.9e9,dene:-9.9e9,ti:-9.9e9,te:-9.9e9,etai:-9.9e9,etae:-9.9e9,b0:-9.9e9}
end

function surf::init,file=file,irho=irho,rhotor=rho,rhomax=rhomax,drhodr=drhodr,rhon=rhon,amin=a,amax=amax,dadr=dadr,an=an,psipol=psi,psimax=psimax,dpsidr=dpsidr,psin=psin,kappa=kappa,dkappadr=dkappadr,delta=delta,ddeltadr=ddeltadr,q=q,dqdr=dqdr,press=p,dpdr=dpdr,d2pdr2=d2pdr2,omega=omega,domegadr=domegadr,rhoi=rhoi,rhoe=rhoe,ni=ni,dene=dene,ti=ti,te=te,etai=etai,etae=etae,b0=b0
;
; CMR August 2003
;
   lfile=keyword_set(file)
   if (not lfile) then begin
      if (keyword_set(irho)) then self.irho=irho
      if (keyword_set(rho)) then self.rho=rho
      if (keyword_set(rhomax)) then self.rhomax=rhomax
      if (keyword_set(drhodr)) then self.drhodr=drhodr
      if (keyword_set(rhon)) then self.rhon=rhon
      if (keyword_set(a)) then self.a=a
      if (keyword_set(amax)) then self.amax=amax
      if (keyword_set(dadr)) then self.dadr=dadr
      if (keyword_set(an)) then self.an=an
      if (keyword_set(psi)) then self.psi=psi
      if (keyword_set(psimax)) then self.psimax=psimax
      if (keyword_set(dpsidr)) then self.dpsidr=dpsidr
      if (keyword_set(psin)) then self.psin=psin
      if (keyword_set(kappa)) then self.kappa=kappa
      if (keyword_set(dkappadr)) then self.dkappadr=dkappadr
      if (keyword_set(delta)) then self.delta=delta
      if (keyword_set(ddeltadr)) then self.ddeltadr=ddeltadr
      if (keyword_set(q)) then self.q=q
      if (keyword_set(dqdr)) then self.dqdr=dqdr
      if (keyword_set(p)) then self.p=p
      if (keyword_set(dpdr)) then self.dpdr=dpdr
      if (keyword_set(d2pdr2)) then self.d2pdr2=d2pdr2
      if (keyword_set(omega)) then self.omega=omega
      if (keyword_set(domegadr)) then self.domegadr=domegadr
      if (keyword_set(rhoi)) then self.rhoi=rhoi
      if (keyword_set(rhoe)) then self.rhoe=rhoe
      if (keyword_set(ni)) then self.ni=ni
      if (keyword_set(dene)) then self.dene=dene
      if (keyword_set(ti)) then self.ti=ti
      if (keyword_set(te)) then self.te=te
      if (keyword_set(etai)) then self.etai=etai
      if (keyword_set(etae)) then self.etae=etae
      if (keyword_set(b0)) then self.b0=b0
   endif else begin
      self->read,file
   endelse
   return,1
end

pro surf::print,lun=lun
;
; CMR August 2003
;
   if (keyword_set(lun)) then begin
      printf,lun,format='("#surfpars"/"#irho",T10,i3,t22,"! flux label for derivatives 1=rhon, 2=an, 3=psin"/"#rho",T10,e12.4,"! m"/"#rhomax",T10,e12.4,"! m"/"#drhodr",T10,e12.4/"#rhon",T10,e12.4/"#a",T10,e12.4,"! m"/"#amax",T10,e12.4,"! m"/"#dadr",T10,e12.4/"#an",T10,e12.4/"#psi",T10,e12.4,"! Wb"/"#psimax",T10,e12.4,"! Wb"/"#dpsidr",T10,e12.4/"#psin",T10,e12.4/"#kappa",T10,e12.4/"#dkappadr",T10,e12.4/"#delta",T10,e12.4/"#ddeltadr",T10,e12.4/"#q",T10,e12.4/"#dqdr",T10,e12.4/"#p",T10,e12.4,"! Pa"/"#dpdr",T10,e12.4/"#d2pdr2",T10,e12.4/,"#omega",T10,e12.4/,"#domegadr",T10,e12.4/"#rhoi",T10,e12.4/"#rhoe",T10,e12.4/"#ni",T10,e12.4/"#dene",T10,e12.4/"#ti",T10,e12.4/"#te",T10,e12.4/"#etai",T10,e12.4/"#etae",T10,e12.4/"#b0",T10,e12.4/"#endsurfpars"/)',self.irho,self.rho,self.rhomax,self.drhodr,self.rhon,self.a,self.amax,self.dadr,self.an,self.psi,self.psimax,self.dpsidr,self.psin,self.kappa,self.dkappadr,self.delta,self.ddeltadr,self.q,self.dqdr,self.p,self.dpdr,self.d2pdr2,self.omega,self.domegadr,self.rhoi,self.rhoe,self.ni,self.dene,self.ti,self.te,self.etai,self.etae,self.b0
   endif else begin
      print,format='("#surfpars"/"#irho",T10,i3,t22,"! flux label for derivatives 1=rhon, 2=an, 3=psin"/"#rho",T10,e12.4,"! m"/"#rhomax",T10,e12.4,"! m"/"#drhodr",T10,e12.4/"#rhon",T10,e12.4/"#a",T10,e12.4,"! m"/"#amax",T10,e12.4,"! m"/"#dadr",T10,e12.4/"#an",T10,e12.4/"#psi",T10,e12.4,"! Wb"/"#psimax",T10,e12.4,"! Wb"/"#dpsidr",T10,e12.4/"#psin",T10,e12.4/"#kappa",T10,e12.4/"#dkappadr",T10,e12.4/"#delta",T10,e12.4/"#ddeltadr",T10,e12.4/"#q",T10,e12.4/"#dqdr",T10,e12.4/"#p",T10,e12.4,"! Pa"/"#dpdr",T10,e12.4/"#d2pdr2",T10,e12.4/,"#omega",T10,e12.4/,"#domegadr",T10,e12.4/"#rhoi",T10,e12.4/"#rhoe",T10,e12.4/"#ni",T10,e12.4/"#dene",T10,e12.4/"#ti",T10,e12.4/"#te",T10,e12.4/"#etai",T10,e12.4/"#etae",T10,e12.4/"#b0",T10,e12.4/"#endsurfpars"/)',self.irho,self.rho,self.rhomax,self.drhodr,self.rhon,self.a,self.amax,self.dadr,self.an,self.psi,self.psimax,self.dpsidr,self.psin,self.kappa,self.dkappadr,self.delta,self.ddeltadr,self.q,self.dqdr,self.p,self.dpdr,self.d2pdr2,self.omega,self.domegadr,self.rhoi,self.rhoe,self.ni,self.dene,self.ti,self.te,self.etai,self.etae,self.b0
   endelse
end

pro surf::read,file
;
; CMR July 2003
;
   get_lun,lun
   a=1.0d0 & i=1L
   openr,lun,file
   rdtill,lun,s,"#surfpars"
   if (strpos(s,'#surfpars') ne -1 ) then begin
      while strpos(s,'#endsurfpars') eq -1 and strpos(s,'#') eq 0 do begin 
         p=strupcase(strtrim(strsplit(s,'# ',/extract),2))
         case p[0] of
         'IRHO': begin
               reads,strmid(s,10),i & self.irho=i
         end
         'RHO': begin
               reads,strmid(s,10),a & self.rho=a
         end
         'RHOMAX': begin
               reads,strmid(s,10),a & self.rhomax=a
         end
         'DRHODR': begin
               reads,strmid(s,10),a & self.drhodr=a
         end
         'RHON': begin
               reads,strmid(s,10),a & self.rhon=a
         end
         'A': begin
               reads,strmid(s,10),a & self.a=a
         end
         'AMAX': begin
               reads,strmid(s,10),a & self.a=amax
         end
         'DADR': begin
               reads,strmid(s,10),a & self.dadr=a
         end
         'AN': begin
               reads,strmid(s,10),a & self.an=a
         end
         'PSI': begin
               reads,strmid(s,10),a & self.psi=a
         end
         'PSIMAX': begin
               reads,strmid(s,10),a & self.psimax=a
         end
         'DPSIDR': begin
               reads,strmid(s,10),a & self.dpsidr=a
         end
         'KAPPA': begin
               reads,strmid(s,10),a & self.kappa=a
         end
         'DKAPPADR': begin
               reads,strmid(s,10),a & self.dkappadr=a
         end
         'DELTA': begin
               reads,strmid(s,10),a & self.delta=a
         end
         'DDELTADR': begin
               reads,strmid(s,10),a & self.ddeltadr=a
         end
         'Q': begin
               reads,strmid(s,10),a & self.q=a
         end
         'DQDR': begin
               reads,strmid(s,10),a & self.dqdr=a
         end
         'P': begin
               reads,strmid(s,10),a & self.p=a
         end
         'DPDR': begin
               reads,strmid(s,10),a & self.dpdr=a
         end
         'D2PDR2': begin
               reads,strmid(s,10),a & self.d2pdr2=a
         end
         'OMEGA': begin
               reads,strmid(s,10),a & self.omega=a
         end
         'DOMEGADR': begin
               reads,strmid(s,10),a & self.domegadr=a
         end
         'RHOI': begin
               reads,strmid(s,10),a & self.rhoi=a
         end
         'RHOE': begin
               reads,strmid(s,10),a & self.rhoe=a
         end
         'NI': begin
               reads,strmid(s,10),a & self.ni=a
         end
         'DENE': begin
               reads,strmid(s,10),a & self.dene=a
         end
         'TI': begin
               reads,strmid(s,10),a & self.ti=a
         end
         'TE': begin
               reads,strmid(s,10),a & self.te=a
         end
         'ETAI': begin
               reads,strmid(s,10),a & self.etai=a
         end
         'ETAE': begin
               reads,strmid(s,10),a & self.etae=a
         end
         'B0': begin
               reads,strmid(s,10),a & self.b0=a
         end
         else:
         endcase
         readf,lun,format='(a80)',s
      endwhile   
   endif
   close,lun
   free_lun,lun
end

function surf::irho
;
; CMR August 2003
;
   return,self.irho
end

function surf::rho
;
; CMR August 2003
;
   return,self.rho
end

function surf::rhomax
;
; CMR August 2003
;
   return,self.rhomax
end

function surf::drhodr
;
; CMR August 2003
;
   return,self.drhodr
end

function surf::rhon
;
; CMR August 2003
;
   return,self.rhon
end

function surf::a
;
; CMR August 2003
;
   return,self.a
end

function surf::dadr
;
; CMR August 2003
;
   return,self.dadr
end

function surf::an
;
; CMR August 2003
;
   return,self.an
end

function surf::psi
;
; CMR August 2003
;
   return,self.psi
end

function surf::psimax
;
; CMR August 2003
;
   return,self.psimax
end

function surf::dpsidr
;
; CMR August 2003
;
   return,self.dpsidr
end

function surf::psin
;
; CMR August 2003
;
   return,self.psin
end

function surf::kappa
;
; CMR August 2003
;
   return,self.kappa
end

function surf::dkappar
;
; CMR August 2003
;
   return,self.dkappar
end

function surf::delta
;
; CMR August 2003
;
   return,self.delta
end

function surf::ddeltar
;
; CMR August 2003
;
   return,self.ddeltar
end

function surf::q
;
; CMR August 2003
;
   return,self.q
end

function surf::dqdr
;
; CMR August 2003
;
   return,self.dqdr
end

function surf::p
;
; CMR August 2003
;
   return,self.p
end

function surf::d2pdr2
;
; CMR August 2003
;
   return,self.d2pdr2
end


function surf::omega
;
; CMR August 2003
;
   return,self.omega
end

function surf::domegadr
;
; CMR August 2003
;
   return,self.domegadr
end

function surf::shearrate
;
; CMR August 2003
;
   pifac=1.6d-19*self.ni*self.ti/self.p
   print,'pifac=',pifac
; suppose irho=3 and all radial derivatives are wrt normalised psi
; assume self.psimax is total poloidal flux (and not per radian)
   if (self.irho ne 3) then begin
      message,string(format='("cannot compute shearrate with irho=",i3)',irho),/continue
   endif
   d2phidr2=self.psimax/(2.0d0*!dpi)*self.domegadr-pifac*self.d2pdr2/(self.ni*1.6d-19)+pifac*self.dpdr^2/(self.ni*1.6d-19*self.p*(1.0d0+self.etai))

   d2phidr2=self.psimax/(2.0d0*!dpi)*self.domegadr-pifac*self.d2pdr2/(self.ni*1.6d-19)+pifac*self.dpdr^2/(self.ni*1.6d-19*self.p*(1.0d0+self.etai))

   wse= 1.0d0/(self.drhodr)^2/self.b0*d2phidr2
   print,'d2phidr2=',d2phidr2
   print,'wse (crude valid only for large aspect ratio)=',wse
   return,wse
end

pro surf::setrhoi,rhoi
;
; CMR August 2003
;
   self.rhoi=rhoi
end

pro surf::setrhoe,rhoe
;
; CMR August 2003
;
   self.rhoe=rhoe
end

pro surf::setni,ni
;
; CMR August 2003
;
   self.ni=ni
end

pro surf::setdene,dene
;
; CMR August 2003
;
   self.dene=dene
end

pro surf::setti,ti
;
; CMR August 2003
;
   self.ti=ti
end

pro surf::sette,te
;
; CMR August 2003
;
   self.te=te
end

pro surf::setetai,etai
;
; CMR August 2003
;
   self.etai=etai
end

pro surf::setetae,etae
;
; CMR August 2003
;
   self.etae=etae
end