mod_thermal_emission.t

Go to the documentation of this file.

! module mod_thermal_emission -- synthesize emission flux of some
! thermal lines
! EUV lines database: 
! 'He_II_304' 'Fe_IX_171' 'Fe_XXIV_193' 'Fe_XIV_211' 'Fe_XVI_335'
! 'Fe_XVIII_94' 'Fe_XXI_131'
! subroutines: 
! get_EUV: get local EUV emission intensity (for 1d, 2d and 3d)
! get_SXR: get local Soft X-ray emission intensity (for 1d, 2d and 3d)
 
module mod_thermal_emission
  use mod_global_parameters
  use mod_geometry
  use mod_physics
  use mod_comm_lib, only: mpistop
 
  implicit none
 
  integer :: n_aia
  double precision :: t_aia(1:101)
  double precision :: f_94(1:101),f_131(1:101),f_171(1:101)
  double precision :: f_193(1:101),f_211(1:101),f_304(1:101)
  double precision :: f_335(1:101)
  integer :: n_iris
  double precision :: t_iris(1:41)
  double precision :: f_1354(1:41)
  integer :: n_eis
  double precision :: t_eis1(1:60),t_eis2(1:60)
  double precision :: f_263(1:60),f_264(1:60),f_192(1:60),f_255(1:60)
 
 
  double precision :: vec_xi1(1:3),vec_xi2(1:3),vec_los(1:3)
 
  data n_aia / 101 /
 
  data t_aia / 4. ,  4.05, 4.1,  4.15, 4.2,  4.25, 4.3,  4.35, &
               4.4,  4.45, 4.5,  4.55, 4.6,  4.65, 4.7,  4.75, &
               4.8,  4.85, 4.9,  4.95, 5. ,  5.05, 5.1,  5.15, &
               5.2,  5.25, 5.3,  5.35, 5.4,  5.45, 5.5,  5.55, &
               5.6,  5.65, 5.7,  5.75, 5.8,  5.85, 5.9,  5.95, &
               6. ,  6.05, 6.1,  6.15, 6.2,  6.25, 6.3,  6.35, &
               6.4,  6.45, 6.5,  6.55, 6.6,  6.65, 6.7,  6.75, &
               6.8,  6.85, 6.9,  6.95, 7. ,  7.05, 7.1,  7.15, &
               7.2,  7.25, 7.3,  7.35, 7.4,  7.45, 7.5,  7.55, &
               7.6,  7.65, 7.7,  7.75, 7.8,  7.85, 7.9,  7.95, &
               8. ,  8.05, 8.1,  8.15, 8.2,  8.25, 8.3,  8.35, &
               8.4,  8.45, 8.5,  8.55, 8.6,  8.65, 8.7,  8.75, &
               8.8,  8.85, 8.9,  8.95, 9. /
 
  data f_94 / 4.25022959d-37, 4.35880298d-36, 3.57054296d-35, 2.18175426d-34, & 
              8.97592571d-34, 2.68512961d-33, 7.49559346d-33, 2.11603751d-32, &
              5.39752853d-32, 1.02935904d-31, 1.33822307d-31, 1.40884290d-31, &
              1.54933156d-31, 2.07543102d-31, 3.42026227d-31, 6.31171444d-31, &
              1.16559416d-30, 1.95360497d-30, 2.77818735d-30, 3.43552578d-30, &
              4.04061803d-30, 4.75470982d-30, 5.65553769d-30, 6.70595782d-30, &
              7.80680354d-30, 8.93247715d-30, 1.02618156d-29, 1.25979030d-29, &
              1.88526483d-29, 3.62448572d-29, 7.50553279d-29, 1.42337571d-28, &
              2.37912813d-28, 3.55232305d-28, 4.84985757d-28, 6.20662827d-28, &
              7.66193687d-28, 9.30403645d-28, 1.10519802d-27, 1.25786927d-27, &
              1.34362634d-27, 1.33185242d-27, 1.22302081d-27, 1.05677973d-27, &
              9.23064720d-28, 8.78570994d-28, 8.02397416d-28, 5.87681142d-28, &
              3.82272695d-28, 3.11492649d-28, 3.85736090d-28, 5.98893519d-28, &
              9.57553548d-28, 1.46650267d-27, 2.10365847d-27, 2.79406671d-27, &
              3.39420087d-27, 3.71077520d-27, 3.57296767d-27, 2.95114380d-27, &
              2.02913103d-27, 1.13361825d-27, 5.13405629d-28, 2.01305089d-28, &
              8.15781482d-29, 4.28366817d-29, 3.08701543d-29, 2.68693906d-29, &
              2.51764203d-29, 2.41773103d-29, 2.33996083d-29, 2.26997246d-29, &
              2.20316143d-29, 2.13810001d-29, 2.07424438d-29, 2.01149189d-29, &
              1.94980213d-29, 1.88917920d-29, 1.82963583d-29, 1.77116920d-29, &
              1.71374392d-29, 1.65740593d-29, 1.60214447d-29, 1.54803205d-29, &
              1.49510777d-29, 1.44346818d-29, 1.39322305d-29, 1.34441897d-29, &
              1.29713709d-29, 1.25132618d-29, 1.20686068d-29, 1.14226584d-29, &
              1.09866413d-29, 1.05635524d-29, 1.01532444d-29, 9.75577134d-30, &
              9.37102736d-30, 8.99873335d-30, 8.63860172d-30, 8.29051944d-30, &
              7.95414793d-30 /
 
  data f_131 / 3.18403601d-37,   3.22254703d-36,   2.61657920d-35, &
               1.59575286d-34,   6.65779556d-34,   2.07015132d-33, &
               6.05768615d-33,   1.76074833d-32,   4.52633001d-32, &
               8.57121883d-32,   1.09184271d-31,   1.10207963d-31, &
               1.11371658d-31,   1.29105226d-31,   1.80385897d-31, &
               3.27295431d-31,   8.92002136d-31,   3.15214579d-30, &
               9.73440787d-30,   2.22709702d-29,   4.01788984d-29, &
               6.27471832d-29,   8.91764995d-29,   1.18725647d-28, &
               1.52888040d-28,   2.05082946d-28,   3.47651873d-28, &
               8.80482184d-28,   2.66533063d-27,   7.05805149d-27, &
               1.46072515d-26,   2.45282476d-26,   3.55303726d-26, &
               4.59075911d-26,   5.36503515d-26,   5.68444094d-26, &
               5.47222296d-26,   4.81119761d-26,   3.85959059d-26, &
               2.80383406d-26,   1.83977650d-26,   1.11182849d-26, &
               6.50748885d-27,   3.96843481d-27,   2.61876319d-27, &
               1.85525324d-27,   1.39717024d-27,   1.11504283d-27, &
               9.38169611d-28,   8.24801234d-28,   7.43331919d-28, &
               6.74537063d-28,   6.14495760d-28,   5.70805277d-28, &
               5.61219786d-28,   6.31981777d-28,   9.19747307d-28, &
               1.76795732d-27,   3.77985446d-27,   7.43166191d-27, &
               1.19785603d-26,   1.48234676d-26,   1.36673114d-26, &
               9.61047146d-27,   5.61209353d-27,   3.04779780d-27, &
               1.69378976d-27,   1.02113491d-27,   6.82223774d-28, &
               5.02099099d-28,   3.99377760d-28,   3.36279037d-28, &
               2.94767378d-28,   2.65740865d-28,   2.44396277d-28, &
               2.28003967d-28,   2.14941419d-28,   2.04178995d-28, &
               1.95031045d-28,   1.87011994d-28,   1.79777869d-28, &
               1.73093957d-28,   1.66795789d-28,   1.60785455d-28, &
               1.55002399d-28,   1.49418229d-28,   1.44022426d-28, &
               1.38807103d-28,   1.33772767d-28,   1.28908404d-28, &
               1.24196208d-28,   1.17437501d-28,   1.12854330d-28, &
               1.08410498d-28,   1.04112003d-28,   9.99529904d-29, &
               9.59358806d-29,   9.20512291d-29,   8.83009123d-29, &
               8.46817043d-29,   8.11921928d-29 /
 
  data f_171 / 2.98015581d-42, 1.24696230d-40, 3.37614652d-39, 5.64103034d-38, &
               5.20550266d-37, 2.77785939d-36, 1.16283616d-35, 6.50007689d-35, &
               9.96177399d-34, 1.89586076d-32, 2.10982799d-31, 1.36946479d-30, &
               6.27396553d-30, 2.29955134d-29, 7.13430211d-29, 1.91024282d-28, &
               4.35358848d-28, 7.94807808d-28, 1.07431875d-27, 1.08399488d-27, &
               9.16212938d-28, 7.34715770d-28, 6.59246382d-28, 9.13541375d-28, &
               2.05939035d-27, 5.08206555d-27, 1.10148083d-26, 2.01884662d-26, &
               3.13578384d-26, 4.14367719d-26, 5.36067711d-26, 8.74170213d-26, &
               1.64161233d-25, 2.94587860d-25, 4.76298332d-25, 6.91765639d-25, &
               9.08825111d-25, 1.08496183d-24, 1.17440114d-24, 1.13943939d-24, &
               9.71696981d-25, 7.09593688d-25, 4.31376399d-25, 2.12708486d-25, &
               8.47429567d-26, 3.17608104d-26, 1.95898842d-26, 1.98064242d-26, &
               1.67706555d-26, 8.99126003d-27, 3.29773878d-27, 1.28896127d-27, &
               8.51169698d-28, 7.53520167d-28, 6.18268143d-28, 4.30034650d-28, &
               2.78152409d-28, 1.95437088d-28, 1.65896278d-28, 1.68740181d-28, &
               1.76054383d-28, 1.63978419d-28, 1.32880591d-28, 1.00833205d-28, &
               7.82252806d-29, 6.36181741d-29, 5.34633869d-29, 4.58013864d-29, &
               3.97833422d-29, 3.49414760d-29, 3.09790940d-29, 2.76786227d-29, &
               2.48806269d-29, 2.24823367d-29, 2.04016653d-29, 1.85977413d-29, &
               1.70367499d-29, 1.56966125d-29, 1.45570643d-29, 1.35964565d-29, &
               1.27879263d-29, 1.21016980d-29, 1.15132499d-29, 1.09959628d-29, &
               1.05307482d-29, 1.01040261d-29, 9.70657096d-30, 9.33214234d-30, &
               8.97689427d-30, 8.63761192d-30, 8.31149879d-30, 7.85162401d-30, &
               7.53828281d-30, 7.23559452d-30, 6.94341530d-30, 6.66137038d-30, &
               6.38929156d-30, 6.12669083d-30, 5.87346434d-30, 5.62943622d-30, & 
               5.39435202d-30 /
 
  data f_193 / 6.40066486d-32, 4.92737300d-31, 2.95342934d-30, 1.28061594d-29, & 
               3.47747667d-29, 5.88554792d-29, 7.72171179d-29, 9.75609282d-29, &
               1.34318963d-28, 1.96252638d-28, 2.70163878d-28, 3.63192965d-28, &
               5.28087341d-28, 8.37821446d-28, 1.39089159d-27, 2.31749718d-27, &
               3.77510689d-27, 5.85198594d-27, 8.26021568d-27, 1.04870405d-26, &
               1.25209374d-26, 1.47406787d-26, 1.77174067d-26, 2.24098537d-26, &
               3.05926105d-26, 4.50018853d-26, 6.84720216d-26, 1.00595861d-25, &
               1.30759222d-25, 1.36481773d-25, 1.15943558d-25, 1.01467304d-25, &
               1.04092532d-25, 1.15071251d-25, 1.27416033d-25, 1.38463476d-25, &
               1.47882726d-25, 1.57041238d-25, 1.69786224d-25, 1.94970397d-25, &
               2.50332918d-25, 3.58321431d-25, 5.18061550d-25, 6.60405549d-25, &
               6.64085365d-25, 4.83825816d-25, 2.40545020d-25, 8.59534098d-26, &
               2.90920638d-26, 1.33204845d-26, 9.03933926d-27, 7.78910836d-27, &
               7.29342321d-27, 7.40267022d-27, 8.05279981d-27, 8.13829291d-27, &
               6.92634262d-27, 5.12521880d-27, 3.59527615d-27, 2.69617560d-27, &
               2.84432713d-27, 5.06697306d-27, 1.01281903d-26, 1.63526978d-26, &
               2.06759342d-26, 2.19482312d-26, 2.10050611d-26, 1.89837248d-26, &
               1.66347131d-26, 1.43071097d-26, 1.21518419d-26, 1.02078343d-26, &
               8.46936184d-27, 6.93015742d-27, 5.56973237d-27, 4.38951754d-27, &
               3.38456457d-27, 2.55309556d-27, 1.88904224d-27, 1.38057546d-27, &
               1.00718330d-27, 7.43581116d-28, 5.63562931d-28, 4.43359435d-28, &
               3.63923535d-28, 3.11248143d-28, 2.75586846d-28, 2.50672237d-28, &
               2.32419348d-28, 2.18325682d-28, 2.06834486d-28, 1.93497044d-28, &
               1.84540751d-28, 1.76356504d-28, 1.68741425d-28, 1.61566157d-28, &
               1.54754523d-28, 1.48249410d-28, 1.42020176d-28, 1.36045230d-28, &
               1.30307965d-28 /
 
  data f_211 / 4.74439912d-42, 1.95251522d-40, 5.19700194d-39, 8.53120166d-38, &
               7.72745727d-37, 4.04158559d-36, 1.64853511d-35, 8.56295439d-35, &
               1.17529722d-33, 2.16867729d-32, 2.40472264d-31, 1.56418133d-30, &
               7.20032889d-30, 2.65838271d-29, 8.33196904d-29, 2.26128236d-28, &
               5.24295811d-28, 9.77791121d-28, 1.35913489d-27, 1.43957785d-27, &
               1.37591544d-27, 1.49029886d-27, 2.06183401d-27, 3.31440622d-27, &
               5.42497318d-27, 8.41100374d-27, 1.17941366d-26, 1.49269794d-26, &
               1.71506074d-26, 1.71266353d-26, 1.51434781d-26, 1.36766622d-26, &
               1.33483562d-26, 1.36834518d-26, 1.45829002d-26, 1.62575306d-26, &
               1.88773347d-26, 2.22026986d-26, 2.54930499d-26, 2.80758138d-26, &
               3.06176409d-26, 3.62799792d-26, 5.13226109d-26, 8.46260744d-26, &
               1.38486586d-25, 1.86192535d-25, 1.78007934d-25, 1.16548409d-25, &
               5.89293257d-26, 2.69952884d-26, 1.24891081d-26, 6.41273176d-27, &
               4.08282914d-27, 3.26463328d-27, 2.76230280d-27, 2.08986882d-27, &
               1.37658470d-27, 8.48489381d-28, 5.19304217d-28, 3.19312514d-28, &
               2.02968197d-28, 1.50171666d-28, 1.39164218d-28, 1.42448821d-28, &
               1.41714519d-28, 1.33341059d-28, 1.20759270d-28, 1.07259692d-28, &
               9.44895400d-29, 8.29030041d-29, 7.25440631d-29, 6.33479483d-29, &
               5.51563757d-29, 4.79002469d-29, 4.14990482d-29, 3.59384972d-29, &
               3.12010860d-29, 2.72624742d-29, 2.40734791d-29, 2.15543565d-29, &
               1.95921688d-29, 1.80682882d-29, 1.68695662d-29, 1.59020936d-29, &
               1.50940886d-29, 1.43956179d-29, 1.37731622d-29, 1.32049043d-29, &
               1.26771875d-29, 1.21803879d-29, 1.17074716d-29, 1.10507836d-29, &
               1.06022834d-29, 1.01703080d-29, 9.75436986d-30, 9.35349257d-30, &
               8.96744546d-30, 8.59527489d-30, 8.23678940d-30, 7.89144480d-30, & 
               7.55891138d-30 /
 
  data f_304 / 3.62695850d-32, 2.79969087d-31, 1.68340584d-30, 7.32681440d-30, &
               1.99967770d-29, 3.41296785d-29, 4.55409104d-29, 5.94994635d-29, &
               8.59864963d-29, 1.39787633d-28, 3.17701965d-28, 1.14474920d-27, &
               4.44845958d-27, 1.54785841d-26, 4.70265345d-26, 1.24524365d-25, &
               2.81535352d-25, 5.10093666d-25, 6.83545307d-25, 6.82110329d-25, &
               5.66886188d-25, 4.36205513d-25, 3.29265688d-25, 2.49802368d-25, &
               1.92527113d-25, 1.51058572d-25, 1.20596047d-25, 9.76884267d-26, &
               7.89979266d-26, 6.18224289d-26, 4.67298332d-26, 3.57934505d-26, &
               2.84535785d-26, 2.32853022d-26, 1.95228514d-26, 1.67880071d-26, &
               1.47608785d-26, 1.32199691d-26, 1.20070960d-26, 1.09378177d-26, &
               1.00031730d-26, 9.62434001d-27, 1.05063954d-26, 1.27267143d-26, &
               1.45923057d-26, 1.36746707d-26, 1.03466970d-26, 6.97647829d-27, &
               4.63141039d-27, 3.19031994d-27, 2.33373613d-27, 1.81589079d-27, &
               1.48446917d-27, 1.26611478d-27, 1.12617468d-27, 1.03625148d-27, &
               9.61400595d-28, 8.79016231d-28, 7.82612130d-28, 6.73762960d-28, &
               5.59717956d-28, 4.53010243d-28, 3.65712196d-28, 3.00958686d-28, &
               2.54011502d-28, 2.18102277d-28, 1.88736437d-28, 1.63817539d-28, &
               1.42283147d-28, 1.23631916d-28, 1.07526003d-28, 9.36797928d-29, &
               8.18565660d-29, 7.18152734d-29, 6.32523238d-29, 5.59513985d-29, &
               4.96614048d-29, 4.42518826d-29, 3.95487628d-29, 3.54690294d-29, &
               3.18953930d-29, 2.87720933d-29, 2.60186750d-29, 2.36011522d-29, &
               2.14717806d-29, 1.95905217d-29, 1.79287981d-29, 1.64562262d-29, &
               1.51489425d-29, 1.39876064d-29, 1.29496850d-29, 1.18665438d-29, &
               1.10240474d-29, 1.02643099d-29, 9.57780996d-30, 8.95465151d-30, &
               8.38950190d-30, 7.87283711d-30, 7.40136507d-30, 6.96804279d-30, & 
               6.56945323d-30 /
 
  data f_335 / 2.46882661d-32, 1.89476632d-31, 1.13216502d-30, 4.89532008d-30, & 
               1.32745970d-29, 2.25390335d-29, 3.00511672d-29, 3.96035934d-29, &
               5.77977656d-29, 8.58600736d-29, 1.14083000d-28, 1.48644411d-28, &
               2.15788823d-28, 3.51628877d-28, 6.12200698d-28, 1.08184987d-27, &
               1.85590697d-27, 2.91679107d-27, 3.94405396d-27, 4.63610680d-27, &
               5.13824456d-27, 5.66602209d-27, 6.30009232d-27, 7.03422868d-27, &
               7.77973918d-27, 8.32371831d-27, 8.56724316d-27, 8.62601374d-27, &
               8.13308844d-27, 6.53188216d-27, 4.55197029d-27, 3.57590087d-27, &
               3.59571707d-27, 4.03502770d-27, 4.54366411d-27, 4.96914990d-27, &
               5.24601170d-27, 5.39979250d-27, 5.43023669d-27, 5.26235042d-27, &
               4.91585495d-27, 4.52628362d-27, 4.13385020d-27, 3.67538967d-27, &
               3.39939742d-27, 3.81284533d-27, 5.02332701d-27, 6.19438602d-27, &
               6.49613071d-27, 6.04010475d-27, 5.24664275d-27, 4.37225997d-27, &
               3.52957182d-27, 2.76212276d-27, 2.08473158d-27, 1.50850518d-27, &
               1.04602472d-27, 7.13091243d-28, 5.34289645d-28, 5.21079581d-28, &
               6.22246365d-28, 6.99555864d-28, 6.29665489d-28, 4.45077026d-28, &
               2.67046793d-28, 1.52774686d-28, 9.18061770d-29, 6.09116074d-29, &
               4.48562572d-29, 3.59463696d-29, 3.05820218d-29, 2.70766652d-29, &
               2.46144034d-29, 2.27758450d-29, 2.13331183d-29, 2.01537836d-29, &
               1.91566180d-29, 1.82893912d-29, 1.75167748d-29, 1.68136168d-29, &
               1.61615595d-29, 1.55481846d-29, 1.49643236d-29, 1.44046656d-29, &
               1.38657085d-29, 1.33459068d-29, 1.28447380d-29, 1.23615682d-29, &
               1.18963296d-29, 1.14478976d-29, 1.10146637d-29, 1.04039479d-29, &
               9.98611410d-30, 9.58205147d-30, 9.19202009d-30, 8.81551313d-30, &
               8.45252127d-30, 8.10224764d-30, 7.76469090d-30, 7.43954323d-30, &
               7.12653873d-30 /
 
 
  data n_iris / 41 /
 
  data t_iris / 4.        , 4.1       , 4.2       , 4.3       , 4.40000001, &
                4.50000001, 4.60000001, 4.70000001, 4.80000001, 4.90000001, &
                5.00000001, 5.10000002, 5.20000002, 5.30000002, 5.40000002, &
                5.50000002, 5.60000002, 5.70000003, 5.80000003, 5.90000003, &
                6.00000003, 6.10000003, 6.20000003, 6.30000003, 6.40000004, &
                6.50000004, 6.60000004, 6.70000004, 6.80000004, 6.90000004, &
                7.00000004, 7.10000005, 7.20000005, 7.30000005, 7.40000005, &
                7.50000005, 7.60000005, 7.70000006, 7.80000006, 7.90000006, &
                8.00000006 /
 
  data f_1354 / 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, &
                0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, &
                0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, &
                0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, &
                0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, &
                0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 1.09503647d-39, &
                5.47214550d-36, 2.42433983d-33, 2.75295034d-31, 1.21929718d-29, &
                2.48392125d-28, 2.33268145d-27, 8.68623633d-27, 1.00166284d-26, &
                3.63126633d-27, 7.45174807d-28, 1.38224064d-28, 2.69270994d-29, &
                5.53314977d-30, 1.15313092d-30, 2.34195788d-31, 4.48242942d-32, &
                7.94976380d-33 /
 
 
  data n_eis  / 60 /
 
  data t_eis1 / 1.99526231d+05, 2.23872114d+05, 2.51188643d+05, 2.81838293d+05, & 
                3.16227766d+05, 3.54813389d+05, 3.98107171d+05, 4.46683592d+05, &
                5.01187234d+05, 5.62341325d+05, 6.30957344d+05, 7.07945784d+05, &
                7.94328235d+05, 8.91250938d+05, 1.00000000d+06, 1.12201845d+06, &
                1.25892541d+06, 1.41253754d+06, 1.58489319d+06, 1.77827941d+06, &
                1.99526231d+06, 2.23872114d+06, 2.51188643d+06, 2.81838293d+06, &
                3.16227766d+06, 3.54813389d+06, 3.98107171d+06, 4.46683592d+06, &
                5.01187234d+06, 5.62341325d+06, 6.30957344d+06, 7.07945784d+06, &
                7.94328235d+06, 8.91250938d+06, 1.00000000d+07, 1.12201845d+07, &
                1.25892541d+07, 1.41253754d+07, 1.58489319d+07, 1.77827941d+07, &
                1.99526231d+07, 2.23872114d+07, 2.51188643d+07, 2.81838293d+07, &
                3.16227766d+07, 3.54813389d+07, 3.98107171d+07, 4.46683592d+07, &
                5.01187234d+07, 5.62341325d+07, 6.30957344d+07, 7.07945784d+07, &
                7.94328235d+07, 8.91250938d+07, 1.00000000d+08, 1.12201845d+08, &
                1.25892541d+08, 1.41253754d+08, 1.58489319d+08, 1.77827941d+08 /
 
  data t_eis2 / 1.99526231d+06, 2.23872114d+06, 2.51188643d+06, 2.81838293d+06, & 
                3.16227766d+06, 3.54813389d+06, 3.98107171d+06, 4.46683592d+06, &
                5.01187234d+06, 5.62341325d+06, 6.30957344d+06, 7.07945784d+06, &
                7.94328235d+06, 8.91250938d+06, 1.00000000d+07, 1.12201845d+07, &
                1.25892541d+07, 1.41253754d+07, 1.58489319d+07, 1.77827941d+07, &
                1.99526231d+07, 2.23872114d+07, 2.51188643d+07, 2.81838293d+07, &
                3.16227766d+07, 3.54813389d+07, 3.98107171d+07, 4.46683592d+07, &
                5.01187234d+07, 5.62341325d+07, 6.30957344d+07, 7.07945784d+07, &
                7.94328235d+07, 8.91250938d+07, 1.00000000d+08, 1.12201845d+08, &
                1.25892541d+08, 1.41253754d+08, 1.58489319d+08, 1.77827941d+08, &
                1.99526231d+08, 2.23872114d+08, 2.51188643d+08, 2.81838293d+08, &
                3.16227766d+08, 3.54813389d+08, 3.98107171d+08, 4.46683592d+08, &
                5.01187234d+08, 5.62341325d+08, 6.30957344d+08, 7.07945784d+08, &
                7.94328235d+08, 8.91250938d+08, 1.00000000d+09, 1.12201845d+09, &
                1.25892541d+09, 1.41253754d+09, 1.58489319d+09, 1.77827941d+09 /
 
  data f_263 /  0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, & 
                0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, &
                0.00000000d+00, 4.46454917d-45, 3.26774829d-42, 1.25292566d-39, &
                2.66922338d-37, 3.28497742d-35, 2.38677554d-33, 1.03937729d-31, &
                2.75075687d-30, 4.47961733d-29, 4.46729177d-28, 2.64862689d-27, &
                8.90863800d-27, 1.72437548d-26, 2.22217752d-26, 2.27999477d-26, &
                2.08264363d-26, 1.78226687d-26, 1.45821699d-26, 1.14675379d-26, &
                8.63082492d-27, 6.15925429d-27, 4.11252514d-27, 2.51530564d-27, &
                1.37090986d-27, 6.42443134d-28, 2.48392636d-28, 7.59187874d-29, &
                1.77852938d-29, 3.23945221d-30, 4.90533903d-31, 6.75458158d-32, &
                9.06878868d-33, 1.23927474d-33, 1.75769395d-34, 2.60710914d-35, &
                4.04318030d-36, 6.53500581d-37, 1.09365022d-37, 1.88383322d-38, &
                3.31425233d-39, 5.90964084d-40, 1.06147549d-40, 1.90706170d-41, &
                3.41331584d-42, 6.07310718d-43, 1.07364738d-43, 1.89085498d-44, &
                3.32598922d-45, 5.87125640d-46, 0.00000000d+00, 0.00000000d+00 /
 
  data f_264 /  0.00000000d+00, 2.81670057d-46, 1.28007268d-43, 2.54586603d-41, & 
                2.67887256d-39, 1.68413285d-37, 6.85702304d-36, 1.91797284d-34, &
                3.84675839d-33, 5.69939170d-32, 6.36224608d-31, 5.39176489d-30, &
                3.45478458d-29, 1.64848693d-28, 5.71476364d-28, 1.39909997d-27, &
                2.37743056d-27, 2.86712530d-27, 2.65206348d-27, 2.07175767d-27, &
                1.47866767d-27, 1.01087374d-27, 6.79605811d-28, 4.54746770d-28, &
                3.04351751d-28, 2.03639149d-28, 1.35940991d-28, 9.01451939d-29, &
                5.91289972d-29, 3.81821178d-29, 2.41434696d-29, 1.48871220d-29, &
                8.93362094d-30, 5.21097445d-30, 2.95964719d-30, 1.64278748d-30, &
                8.95571660d-31, 4.82096011d-31, 2.57390991d-31, 1.36821781d-31, &
                7.27136350d-32, 3.87019426d-32, 2.06883430d-32, 1.11228884d-32, &
                6.01883313d-33, 3.27790676d-33, 1.79805012d-33, 9.93085346d-34, &
                5.52139556d-34, 3.08881387d-34, 1.73890315d-34, 9.84434964d-35, &
                5.60603378d-35, 3.20626492d-35, 1.84111068d-35, 0.00000000d+00, &
                0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00 /
 
  data f_192 /  0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 4.35772105d-44, & 
                1.26162319d-41, 1.97471205d-39, 1.83409019d-37, 1.08206288d-35, &
                4.27914363d-34, 1.17943846d-32, 2.32565755d-31, 3.33087991d-30, &
                3.47013260d-29, 2.60375866d-28, 1.37737127d-27, 5.01053913d-27, &
                1.23479810d-26, 2.11310542d-26, 2.71831513d-26, 2.89851163d-26, &
                2.77312376d-26, 2.50025229d-26, 2.18323661d-26, 1.86980322d-26, &
                1.58035034d-26, 1.31985651d-26, 1.08733133d-26, 8.81804906d-27, &
                7.00417973d-27, 5.43356567d-27, 4.09857884d-27, 2.99651764d-27, &
                2.11902962d-27, 1.45014127d-27, 9.62291023d-28, 6.21548647d-28, &
                3.92807578d-28, 2.44230375d-28, 1.50167782d-28, 9.17611405d-29, &
                5.58707641d-29, 3.40570915d-29, 2.08030862d-29, 1.27588676d-29, &
                7.86535588d-30, 4.87646151d-30, 3.03888897d-30, 1.90578649d-30, &
                1.20195947d-30, 7.61955060d-31, 4.85602199d-31, 3.11049969d-31, &
                2.00087065d-31, 1.29223740d-31, 8.37422008d-32, 0.00000000d+00, &
                0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00 /
 
  data f_255 /  0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 1.76014287d-44, & 
                5.07057938d-42, 7.90473970d-40, 7.31852999d-38, 4.30709255d-36, &
                1.70009061d-34, 4.67925160d-33, 9.21703546d-32, 1.31918676d-30, &
                1.37393161d-29, 1.03102379d-28, 5.45694018d-28, 1.98699648d-27, &
                4.90346776d-27, 8.40524725d-27, 1.08321456d-26, 1.15714525d-26, &
                1.10905152d-26, 1.00155023d-26, 8.75799161d-27, 7.50935839d-27, &
                6.35253533d-27, 5.30919268d-27, 4.37669455d-27, 3.55185164d-27, &
                2.82347055d-27, 2.19257595d-27, 1.65589541d-27, 1.21224987d-27, &
                8.58395132d-28, 5.88163935d-28, 3.90721447d-28, 2.52593407d-28, &
                1.59739995d-28, 9.93802874d-29, 6.11343388d-29, 3.73711135d-29, &
                2.27618743d-29, 1.38793199d-29, 8.48060787d-30, 5.20305940d-30, &
                3.20867365d-30, 1.99011277d-30, 1.24064551d-30, 7.78310544d-31, &
                4.91013681d-31, 3.11338381d-31, 1.98451675d-31, 1.27135460d-31, &
                8.17917486d-32, 5.28280497d-32, 3.42357159d-32, 0.00000000d+00, &
                0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00 /
 
  abstract interface
    subroutine get_subr1(w,x,ixI^L,ixO^L,res)
      use mod_global_parameters
      integer, intent(in)          :: ixI^L, ixO^L
      double precision, intent(in) :: w(ixI^S,nw)
      double precision, intent(in) :: x(ixI^S,1:ndim)
      double precision, intent(out):: res(ixI^S)
    end subroutine get_subr1
 
  end interface
 
  type te_fluid
 
    procedure(get_subr1), pointer, nopass :: get_rho => null()
    procedure(get_subr1), pointer, nopass :: get_pthermal => null()
    procedure(get_subr1), pointer, nopass :: get_var_rfactor => null()
 
  end type te_fluid
 
 
  contains
 
    subroutine get_line_info(wl,ion,mass,logTe,line_center,spatial_px,spectral_px,sigma_PSF,width_slit)
      ! get information of the spectral line
      ! wl: wavelength
      ! mass: ion mass, unit -- proton mass
      ! logTe: peak temperature of emission line in logarithm
      ! line_center: center wavelength of emission line, unit -- Angstrom (0.1 nm) 
      ! spatial_px: pixel size in space of instrument (for image), unit -- arcsec
      ! spectral_px: pixel size in wagelength of instrument (for spectrum), unit -- Angstrom
      ! sigma_PSF: width of point spread function core (for instrument), unit -- pixel
      ! width_slit: width of slit for spectrograph, unit -- arcsec
      use mod_global_parameters
 
      integer, intent(in) :: wl
      integer, intent(out) :: mass
      character(len=30), intent(out) :: ion
      double precision, intent(out) :: logTe,line_center,spatial_px,spectral_px
      double precision, intent(out) :: sigma_PSF,width_slit
 
      select case(wl)
      case(304)
        ion='He II'
        mass=4
        logte=4.7d0
        line_center=303.8d0
        spatial_px=0.6d0
        spectral_px=0.02d0
        sigma_psf=0.895d0
        width_slit=0.6d0
      case(171)
        ion='Fe IX'
        mass=56
        logte=5.8d0
        line_center=171.1d0
        spatial_px=0.6d0
        spectral_px=0.02d0 
        sigma_psf=1.019d0
        width_slit=0.6d0
      case(193)
        ion='Fe XXIV'
        mass=56
        logte=7.3d0
        line_center=193.5d0
        spatial_px=0.6d0
        spectral_px=0.02d0
        sigma_psf=0.813d0
        width_slit=0.6d0
      case(211)
        ion='Fe XIV'
        mass=56
        logte=6.3d0
        line_center=211.3d0
        spatial_px=0.6d0
        spectral_px=0.02d0
        sigma_psf=0.913d0
        width_slit=0.6d0
      case(335)
        ion='Fe XVI'
        mass=56
        logte=6.4d0
        line_center=335.4d0
        spatial_px=0.6d0
        spectral_px=0.02d0
        sigma_psf=1.019d0
        width_slit=0.6d0
      case(94)
        ion='Fe XVIII'
        mass=56
        logte=6.8d0
        line_center=93.9d0
        spatial_px=0.6d0
        spectral_px=0.02d0
        sigma_psf=1.025d0
        width_slit=0.6d0
      case(131)
        ion='Fe XXI'
        mass=56
        logte=7.0d0
        line_center=131.0d0
        spatial_px=0.6d0
        spectral_px=0.02d0
        sigma_psf=0.984d0
        width_slit=0.6d0
      case(1354)
        ion='Fe XXI'
        mass=56
        logte=7.0d0
        line_center=1354.1d0
        spatial_px=0.1663d0
        spectral_px=12.98d-3
        sigma_psf=1.d0
        width_slit=0.33d0
      case(263)
        ion='Fe XVI'
        mass=56
        logte=6.4d0
        line_center=262.976d0
        spatial_px=1.d0
        spectral_px=22.d-3
        sigma_psf=1.d0
        width_slit=2.d0
      case(264)
        ion='Fe XXIII'
        mass=56
        logte=7.1d0
        line_center=263.765d0
        spatial_px=1.d0
        spectral_px=22.d-3
        sigma_psf=1.d0
        width_slit=2.d0
      case(192)
        ion='Fe XXIV'
        mass=56
        logte=7.2d0
        line_center=192.028d0
        spatial_px=1.d0
        spectral_px=22.d-3
        sigma_psf=1.d0
        width_slit=2.d0
      case(255)
        ion='Fe XXIV'
        mass=56
        logte=7.2d0
        line_center=255.113d0
        spatial_px=1.d0
        spectral_px=22.d-3
        sigma_psf=1.d0
        width_slit=2.d0
      case default
        call mpistop("No information about this line")
      end select
 
      spatial_px=spatial_px/instrument_resolution_factor
    end subroutine get_line_info
    
    subroutine get_euv(wl,ixI^L,ixO^L,w,x,fl,flux)
      ! calculate the local emission intensity of given EUV line (optically thin)
      ! wavelength is the wave length of the emission line
      ! unit [DN cm^-1 s^-1 pixel^-1]
      ! ingrate flux along line of sight: DN s^-1 pixel^-1
      use mod_global_parameters
 
      integer, intent(in) :: wl
      integer, intent(in) :: ixI^L, ixO^L
      double precision, intent(in) :: x(ixI^S,1:ndim)
      double precision, intent(in) :: w(ixI^S,1:nw)
      type(te_fluid), intent(in) :: fl
      double precision, intent(out) :: flux(ixI^S)
 
      integer :: n_table
      double precision, allocatable :: t_table(:),f_table(:)
      integer :: ix^D,iTt,i
      double precision :: pth(ixI^S),Te(ixI^S),Ne(ixI^S)
      double precision :: logT,logGT,GT
 
      ! selecting emission table 
      select case(wl)
      case(94)
        n_table=n_aia
        allocate(t_table(1:n_table))
        allocate(f_table(1:n_table))
        t_table(1:n_table)=t_aia(1:n_aia)
        f_table(1:n_table)=f_94(1:n_aia)
      case(131)
        n_table=n_aia
        allocate(t_table(1:n_table))
        allocate(f_table(1:n_table))
        t_table(1:n_table)=t_aia(1:n_aia)
        f_table(1:n_table)=f_131(1:n_aia)
      case(171)
        n_table=n_aia
        allocate(t_table(1:n_table))
        allocate(f_table(1:n_table))
        t_table(1:n_table)=t_aia(1:n_aia)
        f_table(1:n_table)=f_171(1:n_aia)
      case(193)
        n_table=n_aia
        allocate(t_table(1:n_table))
        allocate(f_table(1:n_table))
        t_table(1:n_table)=t_aia(1:n_aia)
        f_table(1:n_table)=f_193(1:n_aia)
      case(211)
        n_table=n_aia
        allocate(t_table(1:n_table))
        allocate(f_table(1:n_table))
        t_table(1:n_table)=t_aia(1:n_aia)
        f_table(1:n_table)=f_211(1:n_aia)
      case(304)
        n_table=n_aia
        allocate(t_table(1:n_table))
        allocate(f_table(1:n_table))
        t_table(1:n_table)=t_aia(1:n_aia)
        f_table(1:n_table)=f_304(1:n_aia)
      case(335)
        n_table=n_aia
        allocate(t_table(1:n_table))
        allocate(f_table(1:n_table))
        t_table(1:n_table)=t_aia(1:n_aia)
        f_table(1:n_table)=f_335(1:n_aia)
      case(1354)
        n_table=n_iris
        allocate(t_table(1:n_table))
        allocate(f_table(1:n_table))
        t_table(1:n_table)=t_iris(1:n_iris)
        f_table(1:n_table)=f_1354(1:n_iris)
      case(263)
        n_table=n_eis
        allocate(t_table(1:n_table))
        allocate(f_table(1:n_table))
        t_table(1:n_table)=t_eis1(1:n_eis)
        f_table(1:n_table)=f_263(1:n_eis)
      case(264)
        n_table=n_eis
        allocate(t_table(1:n_table))
        allocate(f_table(1:n_table))
        t_table(1:n_table)=t_eis2(1:n_eis)
        f_table(1:n_table)=f_264(1:n_eis)
      case(192)
        n_table=n_eis
        allocate(t_table(1:n_table))
        allocate(f_table(1:n_table))
        t_table(1:n_table)=t_eis2(1:n_eis)
        f_table(1:n_table)=f_192(1:n_eis)
      case(255)
        n_table=n_eis
        allocate(t_table(1:n_table))
        allocate(f_table(1:n_table))
        t_table(1:n_table)=t_eis2(1:n_eis)
        f_table(1:n_table)=f_255(1:n_eis)
      case default
        allocate(t_table(1))
        allocate(f_table(1))
        call mpistop("Unknown wavelength")
      end select
      call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
      call fl%get_rho(w,x,ixi^l,ixo^l,ne)
      call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,te)
      te(ixo^s)=pth(ixo^s)/(ne(ixo^s)*te(ixo^s))*unit_temperature
      if (si_unit) then
        ne(ixo^s)=ne(ixo^s)*unit_numberdensity/1.d6 ! m^-3 -> cm-3
        flux(ixo^s)=ne(ixo^s)**2
      else
        ne(ixo^s)=ne(ixo^s)*unit_numberdensity
        flux(ixo^s)=ne(ixo^s)**2
      endif
 
      select case(wl)
      case(94,131,171,193,211,304,335,1354)
      ! temperature table in log
        do i=1,n_table
          if(f_table(i) .gt. 1.d-99) then
            f_table(i)=log10(f_table(i))
          else
            f_table(i)=-99.d0
          endif
        enddo 
        loggt=zero
        {do ix^db=ixomin^db,ixomax^db\}
          logt=log10(te(ix^d))
          if (logt>=t_table(1) .and. logt<=t_table(n_table)) then
            do itt=1,n_table-1
              if (logt>=t_table(itt) .and. logt<t_table(itt+1)) then
                loggt=f_table(itt)*(logt-t_table(itt+1))/(t_table(itt)-t_table(itt+1))+&
                      f_table(itt+1)*(logt-t_table(itt))/(t_table(itt+1)-t_table(itt))
              endif
            enddo
            flux(ix^d)=flux(ix^d)*(10**(loggt))
            if(flux(ix^d)<smalldouble) flux(ix^d)=0.d0
          else
            flux(ix^d)=zero
          endif
        {enddo\}
      case default
      ! temperature table linear
        {do ix^db=ixomin^db,ixomax^db\}
          if (te(ix^d)>=t_table(1) .and. te(ix^d)<=t_table(n_table)) then
            do itt=1,n_table-1
              if (te(ix^d)>=t_table(itt) .and. te(ix^d)<t_table(itt+1)) then
                gt=f_table(itt)*(te(ix^d)-t_table(itt+1))/(t_table(itt)-t_table(itt+1))+&
                   f_table(itt+1)*(te(ix^d)-t_table(itt))/(t_table(itt+1)-t_table(itt))
              endif
            enddo
            flux(ix^d)=flux(ix^d)*gt
            if(flux(ix^d)<smalldouble) flux(ix^d)=0.d0
          else
            flux(ix^d)=zero
          endif
        {enddo\}
      end select
 
      deallocate(t_table,f_table)
    end subroutine get_euv
 
    subroutine get_sxr(ixI^L,ixO^L,w,x,fl,flux,El,Eu)
      !synthesize thermal SXR from El keV to Eu keV released by cm^-3/m^-3
      ! volume of plasma in 1 s
      !flux (cgs): photons cm^-3 s^-1
      !flux (SI): photons m^-3 s^-1
      use mod_global_parameters
 
      integer, intent(in)           :: ixI^L,ixO^L
      integer, intent(in)           :: El,Eu
      double precision, intent(in)  :: x(ixI^S,1:ndim)
      double precision, intent(in)  :: w(ixI^S,nw)
      type(te_fluid), intent(in) :: fl
      double precision, intent(out) :: flux(ixI^S)
 
      integer :: ix^D,ixO^D
      integer :: iE,numE
      double precision :: I0,kb,keV,dE,Ei
      double precision :: pth(ixI^S),Te(ixI^S),kbT(ixI^S)
      double precision :: Ne(ixI^S),gff(ixI^S),fi(ixI^S)
      double precision :: EM(ixI^S)
 
      i0=3.01d-15   ! I0*4*pi*AU**2, I0 from Pinto (2015)
      kb=const_kb
      kev=1.0d3*const_ev
      de=0.1
      nume=floor((eu-el)/de)
      call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
      call fl%get_rho(w,x,ixi^l,ixo^l,ne)
      call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,te)
      te(ixo^s)=pth(ixo^s)/(ne(ixo^s)*te(ixo^s))*unit_temperature
      if (si_unit) then
        ne(ixo^s)=ne(ixo^s)*unit_numberdensity/1.d6 ! m^-3 -> cm-3
        em(ixo^s)=(ne(ixo^s))**2*1.d6 ! cm^-3 m^-3
      else
        ne(ixo^s)=ne(ixo^s)*unit_numberdensity
        em(ixo^s)=(ne(ixo^s))**2
      endif
      kbt(ixo^s)=kb*te(ixo^s)/kev
      flux(ixo^s)=0.0d0
      do ie=0,nume-1
        ei=de*ie+el*1.d0
        gff(ixo^s)=1.d0
        {do ix^db=ixomin^db,ixomax^db\}
          if (kbt(ix^d)>0.01*ei) then
            if(kbt(ix^d)<ei) gff(ix^d)=(kbt(ix^d)/ei)**0.4
            fi(ix^d)=(em(ix^d)*gff(ix^d))*dexp(-ei/(kbt(ix^d)))/(ei*dsqrt(kbt(ix^d)))
          else
            fi(ix^d)=zero
          endif
        {enddo\}
        flux(ixo^s)=flux(ixo^s)+fi(ixo^s)*de
      enddo
      flux(ixo^s)=flux(ixo^s)*i0
    end subroutine get_sxr
 
    subroutine get_goes_sxr_flux(xbox^L,fl,eflux)
      !get GOES SXR 1-8A flux observing at 1AU from given box [w/m^2]
      use mod_global_parameters
 
      double precision, intent(in) :: xbox^L
      type(te_fluid), intent(in) :: fl
      double precision, intent(out) :: eflux
 
      double precision :: dxb^D,xb^L
      integer :: iigrid,igrid,j
      integer :: ixO^L,ixI^L,ix^D
      double precision :: eflux_grid,eflux_pe
 
      ^d&iximin^d=ixglo^d;
      ^d&iximax^d=ixghi^d;
      ^d&ixomin^d=ixmlo^d;
      ^d&ixomax^d=ixmhi^d;
      eflux_pe=zero
      do iigrid=1,igridstail; igrid=igrids(iigrid);
        ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
        ^d&xbmin^d=rnode(rpxmin^d_,igrid);
        ^d&xbmax^d=rnode(rpxmax^d_,igrid);
        call get_goes_flux_grid(ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,ps(igrid)%dvolume(ixi^s),xbox^l,xb^l,fl,eflux_grid)
        eflux_pe=eflux_pe+eflux_grid
      enddo
      call mpi_allreduce(eflux_pe,eflux,1,mpi_double_precision,mpi_sum,icomm,ierrmpi)
    end subroutine get_goes_sxr_flux
 
    subroutine get_goes_flux_grid(ixI^L,ixO^L,w,x,dV,xbox^L,xb^L,fl,eflux_grid)
      use mod_global_parameters
 
      integer, intent(in)           :: ixI^L,ixO^L
      double precision, intent(in)  :: x(ixI^S,1:ndim),dV(ixI^S)
      double precision, intent(in)  :: w(ixI^S,nw)
      double precision, intent(in)  :: xbox^L,xb^L
      type(te_fluid), intent(in) :: fl
      double precision, intent(out) :: eflux_grid
 
      integer :: ix^D,ixO^D,ixb^L
      integer :: iE,numE,j,inbox
      double precision :: I0,kb,keV,dE,Ei,El,Eu,A_cgs
      double precision :: pth(ixI^S),Te(ixI^S),kbT(ixI^S)
      double precision :: Ne(ixI^S),EM(ixI^S)
      double precision :: gff,fi,erg_SI
 
      ! check whether the grid is inside given box
      inbox=0
      {if (xbmin^d<xboxmax^d .and. xbmax^d>xboxmin^d) inbox=inbox+1\}
 
      if (inbox==ndim) then
        ! indexes for cells inside given box
        ^d&ixbmin^d=ixomin^d;
        ^d&ixbmax^d=ixomax^d;
        {if (xbmax^d>xboxmax^d) ixbmax^d=ixomax^d-ceiling((xbmax^d-xboxmax^d)/dxlevel(^d))\}
        {if (xbmin^d<xboxmin^d) ixbmin^d=ceiling((xboxmin^d-xbmin^d)/dxlevel(^d))+ixomin^d\}
 
        i0=1.07d-38 ! photon flux index for observed at 1AU [photon cm^3 m^-2 s^-1 keV^-1]
        kb=const_kb
        kev=1.0d3*const_ev
        erg_si=1.d-7
        a_cgs=1.d-8 ! Angstrom
        el=const_h*const_c/(8.d0*a_cgs)/kev ! 8 A
        eu=const_h*const_c/(1.d0*a_cgs)/kev ! 1 A
        de=0.1  ! keV
        nume=floor((eu-el)/de)
        call fl%get_pthermal(w,x,ixi^l,ixb^l,pth)
        call fl%get_rho(w,x,ixi^l,ixb^l,ne)
        call fl%get_var_Rfactor(w,x,ixi^l,ixb^l,te)
        te(ixb^s)=pth(ixb^s)/(ne(ixb^s)*te(ixb^s))*unit_temperature
        if (si_unit) then
          ne(ixo^s)=ne(ixo^s)*unit_numberdensity/1.d6 ! m^-3 -> cm-3
          em(ixb^s)=(i0*(ne(ixb^s))**2)*dv(ixb^s)*(unit_length*1.d2)**3 ! cm^-3
        else
          ne(ixo^s)=ne(ixo^s)*unit_numberdensity
          em(ixb^s)=(i0*(ne(ixb^s))**2)*dv(ixb^s)*unit_length**3
        endif
        kbt(ixb^s)=kb*te(ixb^s)/kev
        eflux_grid=0.0d0
 
        do ie=0,nume-1
          ei=de*ie+el
          {do ix^db=ixbmin^db,ixbmax^db\}
            if (kbt(ix^d)>1.d-2*ei) then
              if(kbt(ix^d)<ei) then
                gff=(kbt(ix^d)/ei)**0.4
              else
                gff=1.d0
              endif
              fi=(em(ix^d)*gff)*dexp(-ei/(kbt(ix^d)))/(ei*dsqrt(kbt(ix^d)))
              eflux_grid=eflux_grid+fi*de*ei
            endif
          {enddo\}
        enddo
        eflux_grid=eflux_grid*kev*erg_si
      endif
 
    end subroutine get_goes_flux_grid
 
  {^ifthreed
    subroutine get_euv_spectrum(qunit,fl)
      use mod_global_parameters
 
      integer, intent(in) :: qunit
      type(te_fluid), intent(in) :: fl
      character(20) :: datatype
 
      integer :: mass
      character (30) :: ion
      double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
      double precision :: xslit,arcsec
 
      datatype='spectrum_euv'
      arcsec=7.25d7/unit_length
      call get_line_info(spectrum_wl,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
 
      if (mype==0) print *, '###################################################'
      select case(spectrum_wl)
      case (1354)
        if (mype==0) print *, 'Systhesizing EUV spectrum (observed by IRIS).'
      case (263,264,192,255)
        if (mype==0) print *, 'Systhesizing EUV spectrum (observed by Hinode/EIS).'
      case default
        call mpistop('Wrong wavelength!')
      end select
 
      if (spectrum_window_max<=spectrum_window_min) then
        call mpistop('Wrong spectrum window!')
      endif
 
      if (mype==0) write(*,'(a,f8.3,a)') ' Wavelength: ',linecent,' Angstrom'
      if (mype==0) print *, 'Unit of EUV flux: DN s^-1 pixel^-1'
 
      if (dat_resolution) then
        if (mype==0) then
          write(*,'(a,f5.3,a,f5.1,a)') ' Supposed pixel: ',wlrsl,' Angstrom x ',spacersl*725.0, ' km'
          print *, 'Unit of wavelength: Angstrom (0.1 nm) '
          if (si_unit) then
            write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d6,' Mm'
          else
            write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d8,' Mm'
          endif
          write(*,'(a,f8.1,a)') ' Supposed width of slit: ',wslit*725.0,' km'
        endif
        call get_spectrum_datresol(qunit,datatype,fl)
      else
        if (mype==0) then
          print *, 'Unit of wavelength: Angstrom (0.1 nm) '
          if (activate_unit_arcsec) then
            write(*,'(a,f5.3,a,f5.1,a)') ' Pixel: ',wlrsl,' Angstrom x ',spacersl*725.0, ' km'
            print *, 'Unit of length: arcsec (~725 km)'
            write(*,'(a,f8.1,a)') ' Location of slit: xI1 = ',location_slit,' arcsec'
            write(*,'(a,f8.1,a)') ' Width of slit: ',wslit,' arcsec'
          else
            if (si_unit) then
              if (mype==0) write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d6,' Mm'
            else
              if (mype==0) write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d8,' Mm'
            endif
            write(*,'(a,f8.1,a)') ' Location of slit: xI1 = ',location_slit,' Unit_length'
            write(*,'(a,f8.1,a)') ' Width of slit: ',wslit*725.d0,' km'
          endif
        endif
        if (mype==0) print *, 'Direction of the slit: parallel to xI2 vector'
        if (coordinate==cartesian .or. coordinate==spherical) then
          call get_spectrum(qunit,datatype,fl)
        else
          call mpistop("EUV spectrum synthesis: support for sperical coordinates is to be added!")
        endif
      endif
 
      if (mype==0) print *, '###################################################'
 
    end subroutine get_euv_spectrum
 
    subroutine get_spectrum_datresol(qunit,datatype,fl)
 
      integer, intent(in) :: qunit
      character(20), intent(in) :: datatype
      type(te_fluid), intent(in) :: fl
 
      integer :: numWL,numXS,iwL,ixS,numWI,numS
      double precision :: dwLg,xSmin,xSmax,wLmin,wLmax
      double precision, allocatable :: wL(:),xS(:),dwL(:),dxS(:)
      double precision, allocatable :: wI(:,:,:),spectra(:,:),spectra_rc(:,:)
      integer :: strtype,nstrb,nbb,nuni,nstr,bnx
      double precision :: qs,dxfirst,dxmid,lenstr
 
      integer :: iigrid,igrid,j,dir_loc
      double precision :: xbmin(1:ndim),xbmax(1:ndim)
 
      dwlg=1.d-3
      numwl=4*int((spectrum_window_max-spectrum_window_min)/(4.d0*dwlg))
      wlmin=(spectrum_window_max+spectrum_window_min)/2.d0-dwlg*numwl/2
      wlmax=(spectrum_window_max+spectrum_window_min)/2.d0+dwlg*numwl/2
      allocate(wl(numwl),dwl(numwl))
      dwl(:)=dwlg
      do iwl=1,numwl
        wl(iwl)=wlmin+iwl*dwlg-half*dwlg
      enddo
 
      select case(direction_slit)
      case (1)
        numxs=domain_nx1*2**(refine_max_level-1)
        xsmin=xprobmin1
        xsmax=xprobmax1
        bnx=block_nx1
        nbb=domain_nx1
        strtype=stretch_type(1)
        nstrb=nstretchedblocks_baselevel(1)
        qs=qstretch_baselevel(1)
        if (mype==0) print *, 'Direction of the slit: x'
      case (2)
        numxs=domain_nx2*2**(refine_max_level-1)
        xsmin=xprobmin2
        xsmax=xprobmax2
        bnx=block_nx2
        nbb=domain_nx2
        strtype=stretch_type(2)
        nstrb=nstretchedblocks_baselevel(2)
        qs=qstretch_baselevel(2)
        if (mype==0) print *, 'Direction of the slit: y'
      case (3)
        numxs=domain_nx3*2**(refine_max_level-1)
        xsmin=xprobmin3
        xsmax=xprobmax3
        bnx=block_nx3
        nbb=domain_nx3
        strtype=stretch_type(3)
        nstrb=nstretchedblocks_baselevel(3)
        qs=qstretch_baselevel(3)
        if (mype==0) print *, 'Direction of the slit: z'
      case default
        call mpistop('Wrong direction_slit')
      end select
 
      allocate(xs(numxs),dxs(numxs),spectra(numwl,numxs),spectra_rc(numwl,numxs))
      numwi=1
      allocate(wi(numwl,numxs,numwi))
 
      select case(strtype)
      case(0) ! uniform
        dxs(:)=(xsmax-xsmin)/numxs
        do ixs=1,numxs
          xs(ixs)=xsmin+dxs(ixs)*(ixs-half)
        enddo
      case(1) ! uni stretch
        qs=qs**(one/2**(refine_max_level-1))
        dxfirst=(xsmax-xsmin)*(one-qs)/(one-qs**numxs)
        dxs(1)=dxfirst
        do ixs=2,numxs
          dxs(ixs)=dxfirst*qs**(ixs-1)
          xs(ixs)=dxs(1)/(one-qs)*(one-qs**(ixs-1))+half*dxs(ixs)
        enddo
      case(2) ! symm stretch
        ! base level, nbb = nstr + nuni + nstr
        nstr=nstrb*bnx/2
        nuni=nbb-nstrb*bnx
        lenstr=(xsmax-xsmin)/(2.d0+nuni*(one-qs)/(one-qs**nstr))
        dxfirst=(xsmax-xsmin)/(dble(nuni)+2.d0/(one-qs)*(one-qs**nstr))
        dxmid=dxfirst
        ! refine_max level, numXI = nstr + nuni + nstr
        nstr=nstr*2**(refine_max_level-1)
        nuni=nuni*2**(refine_max_level-1)
        qs=qs**(one/2**(refine_max_level-1))
        dxfirst=lenstr*(one-qs)/(one-qs**nstr)
        dxmid=dxmid/2**(refine_max_level-1)
        ! uniform center
        if(nuni .gt. 0) then
          do ixs=nstr+1,nstr+nuni
            dxs(ixs)=dxmid
            xs(ixs)=lenstr+(dble(ixs)-0.5d0-nstr)*dxs(ixs)+xsmin
          enddo
        endif
        ! left half
        do ixs=nstr,1,-1
          dxs(ixs)=dxfirst*qs**(nstr-ixs)
          xs(ixs)=xsmin+lenstr-dxs(ixs)*half-dxfirst*(one-qs**(nstr-ixs))/(one-qs)
        enddo
        ! right half
        do ixs=nstr+nuni+1,numxs
          dxs(ixs)=dxfirst*qs**(ixs-nstr-nuni-1)
          xs(ixs)=xsmax-lenstr+dxs(ixs)*half+dxfirst*(one-qs**(ixs-nstr-nuni-1))/(one-qs)
        enddo
      case default
        call mpistop("unknown stretch type")
      end select
 
      if (los_phi==0 .and. los_theta==90 .and. direction_slit==2) then
      ! LOS->x slit->y
        dir_loc=3
      else if (los_phi==0 .and. los_theta==90 .and. direction_slit==3) then
      ! LOS->x slit->z
        dir_loc=2
      else if (los_phi==90 .and. los_theta==90 .and. direction_slit==1) then
      ! LOS->y slit->x
        dir_loc=3
      else if (los_phi==90 .and. los_theta==90 .and. direction_slit==3) then
      ! LOS->y slit->z
        dir_loc=1
      else if (los_theta==0 .and. direction_slit==1) then
      ! LOS->z slit->x
        dir_loc=2
      else if (los_theta==0 .and. direction_slit==2) then
      ! LOS->z slit->y
        dir_loc=1
      else
        call mpistop('Wrong combination of LOS and slit direction!')
      endif
 
      if (dir_loc==1) then
        if (location_slit>xprobmax1 .or. location_slit<xprobmin1) then
          call mpistop('Wrong value for location_slit!')
        endif
        if(mype==0) write(*,'(a,f8.1,a)') ' Location of slit: x = ',location_slit,' Unit_length'
      else if (dir_loc==2) then
        if (location_slit>xprobmax2 .or. location_slit<xprobmin2) then
          call mpistop('Wrong value for location_slit!')
        endif
        if(mype==0) write(*,'(a,f8.1,a)') ' Location of slit: y = ',location_slit,' Unit_length'
      else
        if (location_slit>xprobmax3 .or. location_slit<xprobmin3) then
          call mpistop('Wrong value for location_slit!')
        endif
        if(mype==0) write(*,'(a,f8.1,a)') ' Location of slit: z = ',location_slit,' Unit_length'
      endif
 
      ! find slit and do integration
      spectra=zero
      do iigrid=1,igridstail; igrid=igrids(iigrid);
        ^d&xbmin(^d)=rnode(rpxmin^d_,igrid);
        ^d&xbmax(^d)=rnode(rpxmax^d_,igrid);
        if (location_slit>=xbmin(dir_loc) .and. location_slit<xbmax(dir_loc)) then
          call integrate_spectra_datresol(igrid,wl,dwl,spectra,numwl,numxs,dir_loc,fl)
        endif
      enddo
 
      nums=numwl*numxs
      call mpi_allreduce(spectra,spectra_rc,nums,mpi_double_precision, &
                         mpi_sum,icomm,ierrmpi)
      do iwl=1,numwl
        do ixs=1,numxs
          if (spectra_rc(iwl,ixs)>smalldouble) then
            wi(iwl,ixs,1)=spectra_rc(iwl,ixs)
          else
            wi(iwl,ixs,1)=zero
          endif
        enddo
      enddo
 
      call output_data(qunit,wl,xs,dwl,dxs,wi,numwl,numxs,numwi,datatype)
 
      deallocate(wl,xs,dwl,dxs,spectra,spectra_rc,wi)
 
    end subroutine get_spectrum_datresol
 
    subroutine integrate_spectra_datresol(igrid,wL,dwL,spectra,numWL,numXS,dir_loc,fl)
      use mod_constants
 
      integer, intent(in) :: igrid,numWL,numXS,dir_loc
      type(te_fluid), intent(in) :: fl
      double precision, intent(in) :: wL(numWL),dwL(numWL)
      double precision, intent(inout) :: spectra(numWL,numXS)
 
      integer :: direction_LOS
      integer :: ixO^L,ixI^L,ix^D,ixOnew
      double precision, allocatable :: flux(:^D&),v(:^D&),pth(:^D&),Te(:^D&),rho(:^D&)
      double precision :: wlc,wlwd
 
      integer :: mass
      double precision :: logTe,lineCent
      character (30) :: ion
      double precision :: spaceRsl,wlRsl,sigma_PSF,wslit
 
      integer :: levelg,rft,ixSmin,ixSmax,iwL
      double precision :: flux_pix,dL
 
      call get_line_info(spectrum_wl,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
 
      if (los_phi==0 .and. los_theta==90) then
        direction_los=1
      else if (los_phi==90 .and. los_theta==90) then
        direction_los=2
      else
        direction_los=3
      endif
 
      ^d&ixomin^d=ixmlo^d\
      ^d&ixomax^d=ixmhi^d\
      ^d&iximin^d=ixglo^d\
      ^d&iximax^d=ixghi^d\
      allocate(flux(ixi^s),v(ixi^s),pth(ixi^s),te(ixi^s),rho(ixi^s))
 
      ^d&ix^d=ixomin^d;
      if (dir_loc==1) then
        do ix1=ixomin1,ixomax1
          if (location_slit>=(ps(igrid)%x(ix^d,1)-half*ps(igrid)%dx(ix^d,1)) .and. &
              location_slit<(ps(igrid)%x(ix^d,1)+half*ps(igrid)%dx(ix^d,1))) then
            ixonew=ix1
          endif
        enddo
        ixomin1=ixonew
        ixomax1=ixonew
      else if (dir_loc==2) then
        do ix2=ixomin2,ixomax2
          if (location_slit>=(ps(igrid)%x(ix^d,2)-half*ps(igrid)%dx(ix^d,2)) .and. &
              location_slit<(ps(igrid)%x(ix^d,2)+half*ps(igrid)%dx(ix^d,2))) then
            ixonew=ix2
          endif
        enddo
        ixomin2=ixonew
        ixomax2=ixonew
      else
        do ix3=ixomin3,ixomax3
          if (location_slit>=(ps(igrid)%x(ix^d,3)-half*ps(igrid)%dx(ix^d,3)) .and. &
              location_slit<(ps(igrid)%x(ix^d,3)+half*ps(igrid)%dx(ix^d,3))) then
            ixonew=ix3
          endif
        enddo
        ixomin3=ixonew
        ixomax3=ixonew
      endif
 
      call get_euv(spectrum_wl,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
      flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2   ! adjust flux due to artifical change of resolution
      call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
      v(ixo^s)=-ps(igrid)%w(ixo^s,iw_mom(direction_los))/rho(ixo^s)
      call fl%get_pthermal(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,pth)
      call fl%get_var_Rfactor(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,te)
      te(ixo^s)=pth(ixo^s)/(te(ixo^s)*rho(ixo^s))
 
      ! grid parameters
      levelg=ps(igrid)%level
      rft=2**(refine_max_level-levelg)
 
      {do ix^d=ixomin^d,ixomax^d\}
        if (flux(ix^d)>smalldouble) then
          if (si_unit) then
            wlc=linecent*(1.d0+v(ix^d)*unit_velocity*1.d2/const_c)
          else
            wlc=linecent*(1.d0+v(ix^d)*unit_velocity/const_c)
          endif
          wlwd=sqrt(kb_cgs*te(ix^d)*unit_temperature/(mass*mp_cgs))
          wlwd=wlwd*linecent/const_c
          ! involved pixel
          select case(direction_slit)
          case(1)
            ixsmin=(block_nx1*(node(pig1_,igrid)-1)+(ix1-ixomin1))*rft+1
            ixsmax=(block_nx1*(node(pig1_,igrid)-1)+(ix1-ixomin1+1))*rft
          case(2)
            ixsmin=(block_nx2*(node(pig2_,igrid)-1)+(ix2-ixomin2))*rft+1
            ixsmax=(block_nx2*(node(pig2_,igrid)-1)+(ix2-ixomin2+1))*rft
          case(3)
            ixsmin=(block_nx3*(node(pig3_,igrid)-1)+(ix3-ixomin3))*rft+1
            ixsmax=(block_nx3*(node(pig3_,igrid)-1)+(ix3-ixomin3+1))*rft
          end select
          ! LOS depth
          select case(direction_los)
          case(1)
            dl=ps(igrid)%dx(ix^d,1)*unit_length
          case(2)
            dl=ps(igrid)%dx(ix^d,2)*unit_length
          case default
            dl=ps(igrid)%dx(ix^d,3)*unit_length
          end select
          if (si_unit) dl=dl*1.d2
          ! integral pixel flux
          do iwl=1,numwl
            flux_pix=flux(ix^d)*wlrsl*dl*exp(-(wl(iwl)-wlc)**2/(2*wlwd**2))/(sqrt(2*dpi)*wlwd)
            if (flux_pix>smalldouble) then
              flux_pix=flux_pix*wslit/spacersl
              spectra(iwl,ixsmin:ixsmax)=spectra(iwl,ixsmin:ixsmax)+flux_pix
            endif
          enddo
        endif
      {enddo\}
 
      deallocate(flux,v,pth,te,rho)
 
    end subroutine integrate_spectra_datresol
 
    subroutine get_spectrum(qunit,datatype,fl)
 
      integer, intent(in) :: qunit
      character(20), intent(in) :: datatype
      type(te_fluid), intent(in) :: fl
 
      integer :: numWL,numXS,iwL,ixS,numWI,ix^D
      double precision :: dwLg,dxSg,xSmin,xSmax,xScent,wLmin,wLmax
      double precision, allocatable :: wL(:),xS(:),dwL(:),dxS(:)
      double precision, allocatable :: wI(:,:,:),spectra(:,:),spectra_rc(:,:)
      double precision :: vec_cor(1:3),xI_cor(1:2)
      double precision :: res,r_loc,r_max
 
      integer :: mass
      character (30) :: ion
      double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
      double precision :: unitv,arcsec,RHESSI_rsl,pixel
      integer :: iigrid,igrid,i,j,numS
      double precision :: xLmin,xLmax,xslit
 
      if (coordinate==spherical) then
        call init_vectors_spherical()
      else
        ! cartesian
        call init_vectors_cartesian()
      endif
 
      ! calculate domain in space
      if (coordinate==spherical) then
        xsmin=-abs(xprobmax1)
        xsmax=abs(xprobmax1)
      else
        do ix1=1,2
          if (ix1==1) vec_cor(1)=xprobmin1
          if (ix1==2) vec_cor(1)=xprobmax1
          do ix2=1,2
            if (ix2==1) vec_cor(2)=xprobmin2
            if (ix2==2) vec_cor(2)=xprobmax2
            do ix3=1,2
              if (ix3==1) vec_cor(3)=xprobmin3
              if (ix3==2) vec_cor(3)=xprobmax3
              if (big_image) then
                r_loc=(vec_cor(1)-x_origin(1))**2
                r_loc=r_loc+(vec_cor(2)-x_origin(2))**2
                r_loc=r_loc+(vec_cor(3)-x_origin(3))**2
                r_loc=sqrt(r_loc)
                if (ix1==1 .and. ix2==1 .and. ix3==1) then
                  r_max=r_loc
                else
                  r_max=max(r_max,r_loc)
                endif
              else
                call get_cor_image(vec_cor,xi_cor)
                if (ix1==1 .and. ix2==1 .and. ix3==1) then
                  xsmin=xi_cor(2)
                  xsmax=xi_cor(2)
                else
                  xsmin=min(xsmin,xi_cor(2))
                  xsmax=max(xsmax,xi_cor(2))
                endif
              endif
            enddo
          enddo
        enddo
        if (big_image) then
          xsmin=-r_max
          xsmax=r_max
        endif
      endif
      xscent=(xsmin+xsmax)/2.d0
 
      ! tables for storing spectra data
      if (si_unit) then
        arcsec=7.25d5/unit_length
      else
        arcsec=7.25d7/unit_length
      endif
      call get_line_info(spectrum_wl,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
      dxsg=spacersl*arcsec
      numxs=ceiling((xsmax-xscent)/dxsg)
      xsmin=xscent-numxs*dxsg
      xsmax=xscent+numxs*dxsg
      numxs=numxs*2
      dwlg=wlrsl
      numwl=2*int((spectrum_window_max-spectrum_window_min)/(2.d0*dwlg))
      wlmin=(spectrum_window_max+spectrum_window_min)/2.d0-dwlg*numwl/2
      wlmax=(spectrum_window_max+spectrum_window_min)/2.d0+dwlg*numwl/2
      allocate(wl(numwl),dwl(numwl),xs(numxs),dxs(numxs))
      numwi=1
      allocate(wi(numwl,numxs,numwi),spectra(numwl,numxs),spectra_rc(numwl,numxs))
      do iwl=1,numwl
        wl(iwl)=wlmin+iwl*dwlg-half*dwlg
        dwl=dwlg
      enddo
      do ixs=1,numxs
        xs(ixs)=xsmin+dxsg*(ixs-half)
        dxs(ixs)=dxsg
      enddo
 
      ! find slit and do integration
      spectra=zero
      do iigrid=1,igridstail; igrid=igrids(iigrid);
        do ix1=1,2
          if (ix1==1) vec_cor(1)=rnode(rpxmin1_,igrid)
          if (ix1==2) vec_cor(1)=rnode(rpxmax1_,igrid)
          do ix2=1,2
            if (ix2==1) vec_cor(2)=rnode(rpxmin2_,igrid)
            if (ix2==2) vec_cor(2)=rnode(rpxmax2_,igrid)
            do ix3=1,2
              if (ix3==1) vec_cor(3)=rnode(rpxmin3_,igrid)
              if (ix3==2) vec_cor(3)=rnode(rpxmax3_,igrid)
              call get_cor_image(vec_cor,xi_cor)
              if (ix1==1 .and. ix2==1 .and. ix3==1) then
                xlmin=xi_cor(1)
                xlmax=xi_cor(1)
              else
                xlmin=min(xlmin,xi_cor(1))
                xlmax=max(xlmax,xi_cor(1))
              endif
            enddo
          enddo
        enddo
 
        if (activate_unit_arcsec) then 
          xslit=location_slit*arcsec
        else
          xslit=location_slit
        endif
        if (xslit>=xlmin-wslit*arcsec .and. xslit<=xlmax+wslit*arcsec) then
          call integrate_spectra_cartesian(igrid,wl,dwlg,xs,dxsg,spectra,numwl,numxs,fl)
        endif
      enddo
 
      nums=numwl*numxs
      call mpi_allreduce(spectra,spectra_rc,nums,mpi_double_precision, &
                         mpi_sum,icomm,ierrmpi)
      do iwl=1,numwl
        do ixs=1,numxs
          if (spectra_rc(iwl,ixs)>smalldouble) then
            wi(iwl,ixs,1)=spectra_rc(iwl,ixs)
          else
            wi(iwl,ixs,1)=zero
          endif
        enddo
      enddo
 
      if (activate_unit_arcsec) then 
        xs=xs/arcsec
        dxs=dxs/arcsec
      endif
 
      call output_data(qunit,wl,xs,dwl,dxs,wi,numwl,numxs,numwi,datatype)
 
      deallocate(wl,xs,dwl,dxs,spectra,spectra_rc,wi)
 
    end subroutine get_spectrum
 
    subroutine integrate_spectra_cartesian(igrid,wL,dwLg,xS,dxSg,spectra,numWL,numXS,fl)
 
      integer, intent(in) :: igrid,numWL,numXS
      double precision, intent(in) :: wL(numWL),xS(numXS)
      double precision, intent(in) :: dwLg,dxSg
      double precision, intent(inout) :: spectra(numWL,numXS)
      type(te_fluid), intent(in) :: fl
 
      integer :: ixO^L,ixI^L,ix^D,ixOnew,j
      double precision, allocatable :: flux(:^D&),v(:^D&),pth(:^D&),Te(:^D&),rho(:^D&)
      double precision :: wlc,wlwd,res,dst_slit,xslit,arcsec
      double precision :: vloc(1:3),xloc(1:3),dxloc(1:3),xIloc(1:2),dxIloc(1:2)
      integer :: nSubC^D,iSubC^D,iwL,ixS,ixSmin,ixSmax,iwLmin,iwLmax,nwL
      double precision :: slit_width,dxSubC^D,xerf^L,fluxSubC
      double precision :: xSubC(1:3),xCent(1:2)
 
      integer :: mass
      double precision :: logTe,lineCent
      character (30) :: ion
      double precision :: spaceRsl,wlRsl,sigma_PSF,wslit
      double precision :: sigma_wl,sigma_xs,factor
 
      if (si_unit) then
        arcsec=7.25d5/unit_length
      else
        arcsec=7.25d7/unit_length
      endif
      if (activate_unit_arcsec) then 
        xslit=location_slit*arcsec
      else
        xslit=location_slit
      endif
 
      call get_line_info(spectrum_wl,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
 
      ^d&ixomin^d=ixmlo^d\
      ^d&ixomax^d=ixmhi^d\
      ^d&iximin^d=ixglo^d\
      ^d&iximax^d=ixghi^d\
      allocate(flux(ixi^s),v(ixi^s),pth(ixi^s),te(ixi^s),rho(ixi^s))
      ! get local EUV flux and velocity
      call get_euv(spectrum_wl,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
      flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2   ! adjust flux due to artifical change of resolution
      call fl%get_pthermal(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,pth)
      call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
      call fl%get_var_Rfactor(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,te)
      te(ixo^s)=pth(ixo^s)/(te(ixo^s)*rho(ixo^s))
      {do ix^d=ixomin^d,ixomax^d\}
        do j=1,3
          vloc(j)=ps(igrid)%w(ix^d,iw_mom(j))/rho(ix^d)
        enddo
        call dot_product_loc(vloc,vec_los,res)
        v(ix^d)=res
      {enddo\}
 
      deallocate(rho)
 
      slit_width=wslit*arcsec
      sigma_wl=sigma_psf*dwlg
      sigma_xs=sigma_psf*dxsg
      {do ix^d=ixomin^d,ixomax^d\}
        if (flux(ix^d)>smalldouble) then
          xloc(1:3)=ps(igrid)%x(ix^d,1:3)
          dxloc(1:3)=ps(igrid)%dx(ix^d,1:3)
          call get_cor_image(xloc,xiloc)
          call dot_product_loc(dxloc,vec_xi1,res)
          dxiloc(1)=abs(res)
          if (xiloc(1)>=xslit-half*(slit_width+dxiloc(1)) .and. & 
              xiloc(1)<=xslit+half*(slit_width+dxiloc(1))) then
            ^d&nsubc^d=1;
            ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(vec_xi1(^d))/(slit_width/16.d0)));
            ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(vec_xi2(^d))/(dxsg/4.d0)));
            ^d&dxsubc^d=ps(igrid)%dx(ix^dd,^d)/nsubc^d;
            ! local line center and line width
            if (si_unit) then
              fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length*1.d2/dxsg/dxsg  ! DN s^-1
              wlc=linecent*(1.d0+v(ix^d)*unit_velocity*1.d2/const_c)
            else
              fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length/dxsg/dxsg  ! DN s^-1
              wlc=linecent*(1.d0+v(ix^d)*unit_velocity/const_c)
            endif
            wlwd=sqrt(kb_cgs*te(ix^d)*unit_temperature/(mass*mp_cgs))
            wlwd=wlwd*linecent/const_c
            ! dividing a cell to several parts to get more accurate integrating values
            {do isubc^d=1,nsubc^d\}
              ^d&xsubc(^d)=xloc(^d)-half*dxloc(^d)+(isubc^d-half)*dxsubc^d;
              call get_cor_image(xsubc,xcent)
              dst_slit=abs(xcent(1)-xslit)  ! space distance to slit center
              if (dst_slit<=half*slit_width) then
                ixs=floor((xcent(2)-(xs(1)-half*dxsg))/dxsg)+1
                ixsmin=max(1,ixs-3)
                ixsmax=min(ixs+3,numxs)
                iwl=floor((wlc-(wl(1)-half*dwlg))/dwlg)+1
                nwl=3*ceiling(wlwd/dwlg+1)
                iwlmin=max(1,iwl-nwl)
                iwlmax=min(iwl+nwl,numwl)
                ! calculate the contribution to nearby pixels
                do iwl=iwlmin,iwlmax
                  do ixs=ixsmin,ixsmax
                    xerfmin1=(wl(iwl)-half*dwlg-wlc)/sqrt(2.d0*(sigma_wl**2+wlwd**2))
                    xerfmax1=(wl(iwl)+half*dwlg-wlc)/sqrt(2.d0*(sigma_wl**2+wlwd**2))
                    xerfmin2=(xs(ixs)-half*dxsg-xcent(2))/(sqrt(2.d0)*sigma_xs)
                    xerfmax2=(xs(ixs)+half*dxsg-xcent(2))/(sqrt(2.d0)*sigma_xs)
                    factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
                    spectra(iwl,ixs)=spectra(iwl,ixs)+fluxsubc*factor
                  enddo
                enddo
                ! nearby pixels
              endif
            {enddo\}
          endif
        endif
      {enddo\}
 
      deallocate(flux,v,pth,te)
    end subroutine integrate_spectra_cartesian
  }
 
  {^ifthreed
    subroutine get_euv_image(qunit,fl)
      use mod_global_parameters
 
      integer, intent(in) :: qunit
      type(te_fluid), intent(in) :: fl
      character(20) :: datatype
 
      integer :: mass
      character (30) :: ion
      double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
      double precision :: t0,t1
 
      t0=mpi_wtime()
      datatype='image_euv'
      call get_line_info(wavelength,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
 
      if (mype==0) then
        print *, '###################################################'
        print *, 'Systhesizing EUV image'
        write(*,'(a,f8.3,a)') ' Wavelength: ',linecent,' Angstrom'
        print *, 'Unit of EUV flux: DN s^-1 pixel^-1'
      endif
 
      if (dat_resolution) then
        if (coordinate/=cartesian) call mpistop('EUV synthesis: only cartesian is supported for .dat resolution!')
        if (mype==0) then
          write(*,'(a,f7.1,a,f7.1,a,f5.1,a,f5.1,a)') ' Supposed Pixel: ',spacersl*725.0,' km x ',spacersl*725.0, & 
                                                     ' km  (', spacersl, ' arcsec x ', spacersl, ' arcsec)'
          if (si_unit) then
            write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d6,' Mm'
          else
            write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d8,' Mm'
          endif
        endif
        if (los_phi==0 .and. los_theta==90) then
          call get_image_datresol(qunit,datatype,fl)
        else if (los_phi==90 .and. los_theta==90) then
          call get_image_datresol(qunit,datatype,fl)
        else if (los_theta==0) then
          call get_image_datresol(qunit,datatype,fl)
        else
          call mpistop('ERROR: Wrong LOS for synthesizing emission!')
        endif
      else
        if (mype==0) then
          write(*,'(a,f7.1,a,f7.1,a,f5.1,a,f5.1,a)') ' Pixel: ',spacersl*725.0,' km x ',spacersl*725.0, ' km  (', &
                                                                  spacersl, ' arcsec x ', spacersl, ' arcsec)'
          if (activate_unit_arcsec) then
            print *, 'Unit of length: arcsec (~725 km)'
          else
            if (si_unit) then
              write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d6,' Mm'
            else
              write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d8,' Mm'
            endif
          endif
        endif
        if (coordinate==cartesian) then
          if (mype==0) write(*,'(a,f8.3,f8.3,f8.3,a)') ' Mapping: [',x_origin(1),x_origin(2),x_origin(3), &
                                                     '] of the simulation box is located at [X=0,Y=0] of the image'
            call get_image(qunit,datatype,fl)
        else if (coordinate==spherical) then
          if (mype==0) write(*,'(a,f6.3,f8.3,f8.3,a)') ' Mapping: R=0 of the simulation box is located at [X=0,Y=0] of the image'
          call get_image(qunit,datatype,fl)
        else
          call mpistop("EUV synthesis: this coordinate is not supported!")
        endif
      endif
 
      t1=mpi_wtime()
      if (mype==0) print *, 'time comsuming: ',t1-t0,' s'
      if (mype==0) print *, '###################################################'
      
    end subroutine get_euv_image
 
    subroutine get_sxr_image(qunit,fl)
      use mod_global_parameters
 
      integer, intent(in) :: qunit
      type(te_fluid), intent(in) :: fl
      character(20) :: datatype
      double precision :: RHESSI_rsl
      double precision :: t0,t1
 
      t0=mpi_wtime()
      datatype='image_sxr'
      rhessi_rsl=2.3/instrument_resolution_factor
 
      if (mype==0) then
        print *, '###################################################'
        print *, 'Systhesizing SXR image (observed at 1 AU).'
        write(*,'(a,i2,a,i2,a)') ' Passband: ',emin_sxr,' - ',emax_sxr,' keV'
      endif
 
      if (dat_resolution) then
        if (coordinate/=cartesian) call mpistop('SXR synthesis: only cartesian is supported for .dat resolution!')
        if (mype==0) then
          print *, 'Unit of SXR flux: photons cm^-2 s^-1 pixel^-1'
          write(*,'(a,f5.1,a,f5.1,a,f5.1,a,f5.1,a)') ' Supposed Pixel: ',rhessi_rsl*0.725, ' Mm x ',rhessi_rsl*0.725, &
                                                     ' Mm  (', rhessi_rsl, ' arcsec x ', rhessi_rsl, ' arcsec)'
          if (si_unit) then
            write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d6,' Mm'
          else
            write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d8,' Mm'
          endif
        endif
        if (los_phi==0 .and. los_theta==90) then
          call get_image_datresol(qunit,datatype,fl)
        else if (los_phi==90 .and. los_theta==90) then
          call get_image_datresol(qunit,datatype,fl)
        else if (los_theta==0) then
          call get_image_datresol(qunit,datatype,fl)
        else
          call mpistop('ERROR: Wrong LOS for synthesizing emission!')
        endif
      else
        if (mype==0) then
          print *, 'Unit of SXR flux: photons cm^-2 s^-1 pixel^-1'
          write(*,'(a,f5.1,a,f5.1,a,f5.1,a,f5.1,a)') ' Pixel: ',rhessi_rsl*0.725, ' Mm x ',rhessi_rsl*0.725, &
                                                      ' Mm  (', rhessi_rsl, ' arcsec x ', rhessi_rsl, ' arcsec)'
          if (activate_unit_arcsec) then
            print *, 'Unit of length: arcsec (~725 km)'
          else
            if (si_unit) then
              write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d6,' Mm'
            else
              write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d8,' Mm'
            endif
          endif
        endif
        if (coordinate==cartesian) then
          if (mype==0) write(*,'(a,f8.3,f8.3,f8.3,a)') ' Mapping: [',x_origin(1),x_origin(2),x_origin(3), &
                                                     '] of the simulation box is located at [X=0,Y=0] of the image'
            call get_image(qunit,datatype,fl)
        else if (coordinate==spherical) then
          if (mype==0) write(*,'(a,f6.3,f8.3,f8.3,a)') ' Mapping: R=0 of the simulation box is located at [X=0,Y=0] of the image'
          call get_image(qunit,datatype,fl)
        else
          call mpistop("SXR synthesis: this coordinate is not supported!")
        endif
      endif
 
      t1=mpi_wtime()
      if (mype==0) print *, 'time comsuming:',t1-t0
      if (mype==0) print *, '###################################################'
 
    end subroutine get_sxr_image
 
    subroutine get_whitelight_image(qunit,fl)
      use mod_global_parameters
 
      integer, intent(in) :: qunit
      type(te_fluid), intent(in) :: fl
      character(20) :: datatype
      double precision :: LASCO_rsl
 
      if (mype==0) print *, '###################################################'
 
      if (whitelight_instrument=='LASCO/C1') then
        lasco_rsl=5.6d0/instrument_resolution_factor
        if (mype==0) print *, 'Systhesizing white light image (observed by LASCO/C1).'
      else if (whitelight_instrument=='LASCO/C2') then
        lasco_rsl=11.4d0/instrument_resolution_factor
        if (mype==0) print *, 'Systhesizing white light image (observed by LASCO/C2).'
      else if (whitelight_instrument=='LASCO/C3') then
        lasco_rsl=56.d0/instrument_resolution_factor
        if (mype==0) print *, 'Systhesizing white light image (observed by LASCO/C3).'
      else
        call mpistop('Whitelight synthesis: instrument is not supported!')
      endif
 
      if (mype==0) write(*,'(a,f5.1,a,f5.1,a,f5.1,a,f5.1,a)') ' Pixel: ',lasco_rsl*0.725,' Mm x ',lasco_rsl*0.725, ' Mm  (', &
                                                                lasco_rsl, ' arcsec x ', lasco_rsl, ' arcsec) '
      if (mype==0) print *, 'Unit of white light flux: average Sun brightness'
 
      datatype='image_whitelight'
 
      if (mype==0) then
        if (activate_unit_arcsec) then
          print *, 'Unit of length: arcsec (~725 km)'
        else
          if (si_unit) then
            if (mype==0) write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d6,' Mm'
          else
            if (mype==0) write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d8,' Mm'
          endif
        endif
      endif
 
      if (coordinate==spherical) then
        if (mype==0) write(*,'(a,f6.3,f8.3,f8.3,a)') ' Mapping: R=0 of the simulation box is located at [X=0,Y=0] of the image'
        call get_image(qunit,datatype,fl)
      else
        call mpistop("Whitelight synthesis: this coordinate is not supported!")
      endif
 
      if (mype==0) print *, '###################################################'
 
    end subroutine get_whitelight_image
 
    subroutine get_image_datresol(qunit,datatype,fl)
      ! integrate emission flux along line of sight (LOS) 
      ! in a 3D simulation box and get a 2D EUV image
      use mod_global_parameters
      use mod_constants
 
      integer, intent(in) :: qunit
      character(20), intent(in) :: datatype
      type(te_fluid), intent(in) :: fl
 
      double precision :: dx^D
      integer :: numX^D,ix^D
      double precision, allocatable :: EUV(:,:),EUVs(:,:),Dpl(:,:),Dpls(:,:)
      double precision, allocatable :: SXR(:,:),SXRs(:,:),wI(:,:,:)
      double precision, allocatable :: xI1(:),xI2(:),dxI1(:),dxI2(:),dxIi
      integer :: numXI1,numXI2,numSI,numWI
      double precision :: xI^L
      integer :: iigrid,igrid,i,j
      double precision, allocatable :: xIF1(:),xIF2(:),dxIF1(:),dxIF2(:)
      integer :: nXIF1,nXIF2
      double precision :: xIF^L
 
      double precision :: unitv,arcsec,RHESSI_rsl
      integer :: strtype^D,nstrb^D,nbb^D,nuni^D,nstr^D,bnx^D
      double precision :: qs^D,dxfirst^D,dxmid^D,lenstr^D
 
      numx1=domain_nx1*2**(refine_max_level-1)
      numx2=domain_nx2*2**(refine_max_level-1)
      numx3=domain_nx3*2**(refine_max_level-1)
 
      ! parameters for creating table
      if (los_phi==0 .and. los_theta==90) then
        nxif1=domain_nx2*2**(refine_max_level-1)
        nxif2=domain_nx3*2**(refine_max_level-1)
        xifmin1=xprobmin2
        xifmax1=xprobmax2
        xifmin2=xprobmin3
        xifmax2=xprobmax3
        bnx1=block_nx2
        bnx2=block_nx3
        nbb1=domain_nx2
        nbb2=domain_nx3
        strtype1=stretch_type(2)
        strtype2=stretch_type(3)
        nstrb1=nstretchedblocks_baselevel(2)
        nstrb2=nstretchedblocks_baselevel(3)
        qs1=qstretch_baselevel(2)
        qs2=qstretch_baselevel(3)
        if (mype==0) write(*,'(a)') ' LOS vector: [-1.00  0.00  0.00]'
        if (mype==0) write(*,'(a)') ' xI1 vector: [ 0.00  1.00  0.00]'
        if (mype==0) write(*,'(a)') ' xI2 vector: [ 0.00  0.00  1.00]'
      else if (los_phi==90 .and. los_theta==90) then
        nxif1=domain_nx3*2**(refine_max_level-1)
        nxif2=domain_nx1*2**(refine_max_level-1)
        xifmin1=xprobmin3
        xifmax1=xprobmax3
        xifmin2=xprobmin1
        xifmax2=xprobmax1
        bnx1=block_nx3
        bnx2=block_nx1
        nbb1=domain_nx3
        nbb2=domain_nx1
        strtype1=stretch_type(3)
        strtype2=stretch_type(1)
        nstrb1=nstretchedblocks_baselevel(3)
        nstrb2=nstretchedblocks_baselevel(1)
        qs1=qstretch_baselevel(3)
        qs2=qstretch_baselevel(1)
        if (mype==0) write(*,'(a)') ' LOS vector: [ 0.00 -1.00  0.00]'
        if (mype==0) write(*,'(a)') ' xI1 vector: [-1.00  0.00  0.00]'
        if (mype==0) write(*,'(a)') ' xI2 vector: [ 0.00  0.00  1.00]'
      else
        nxif1=domain_nx1*2**(refine_max_level-1)
        nxif2=domain_nx2*2**(refine_max_level-1)
        xifmin1=xprobmin1
        xifmax1=xprobmax1
        xifmin2=xprobmin2
        xifmax2=xprobmax2
        bnx1=block_nx1
        bnx2=block_nx2
        nbb1=domain_nx1
        nbb2=domain_nx2
        strtype1=stretch_type(1)
        strtype2=stretch_type(2)
        nstrb1=nstretchedblocks_baselevel(1)
        nstrb2=nstretchedblocks_baselevel(2)
        qs1=qstretch_baselevel(1)
        qs2=qstretch_baselevel(2)
        if (mype==0) write(*,'(a)') ' LOS vector: [ 0.00  0.00 -1.00]'
        if (mype==0) write(*,'(a)') ' xI1 vector: [ 1.00  0.00  0.00]'
        if (mype==0) write(*,'(a)') ' xI2 vector: [ 0.00  1.00  0.00]'
      endif
      allocate(xif1(nxif1),xif2(nxif2),dxif1(nxif1),dxif2(nxif2))
 
      ! initialize image coordinate
      select case(strtype1)
      case(0) ! uniform
        dxif1(:)=(xifmax1-xifmin1)/nxif1
        do ix1=1,nxif1
          xif1(ix1)=xifmin1+dxif1(ix1)*(ix1-half)
        enddo
      case(1) ! uni stretch
        qs1=qs1**(one/2**(refine_max_level-1))
        dxfirst1=(xifmax1-xifmin1)*(one-qs1)/(one-qs1**nxif1)
        dxif1(1)=dxfirst1
        do ix1=2,nxif1
          dxif1(ix1)=dxfirst1*qs1**(ix1-1)
          xif1(ix1)=dxif1(1)/(one-qs1)*(one-qs1**(ix1-1))+half*dxif1(ix1)
        enddo
      case(2) ! symm stretch
        ! base level, nbb = nstr + nuni + nstr
        nstr1=nstrb1*bnx1/2
        nuni1=nbb1-nstrb1*bnx1
        lenstr1=(xifmax1-xifmin1)/(2.d0+nuni1*(one-qs1)/(one-qs1**nstr1))
        dxfirst1=(xifmax1-xifmin1)/(dble(nuni1)+2.d0/(one-qs1)*(one-qs1**nstr1))
        dxmid1=dxfirst1
        ! refine_max level, numXI = nstr + nuni + nstr
        nstr1=nstr1*2**(refine_max_level-1)
        nuni1=nuni1*2**(refine_max_level-1)
        qs1=qs1**(one/2**(refine_max_level-1))
        dxfirst1=lenstr1*(one-qs1)/(one-qs1**nstr1)
        dxmid1=dxmid1/2**(refine_max_level-1)
        ! uniform center
        if(nuni1 .gt. 0) then
          do ix1=nstr1+1,nstr1+nuni1
            dxif1(ix1)=dxmid1
            xif1(ix1)=lenstr1+(dble(ix1)-0.5d0-nstr1)*dxif1(ix1)+xifmin1
          enddo
        endif
        ! left half
        do ix1=nstr1,1,-1
          dxif1(ix1)=dxfirst1*qs1**(nstr1-ix1)
          xif1(ix1)=xifmin1+lenstr1-dxif1(ix1)*half-dxfirst1*(one-qs1**(nstr1-ix1))/(one-qs1)
        enddo
        ! right half
        do ix1=nstr1+nuni1+1,nxif1
          dxif1(ix1)=dxfirst1*qs1**(ix1-nstr1-nuni1-1)
          xif1(ix1)=xifmax1-lenstr1+dxif1(ix1)*half+dxfirst1*(one-qs1**(ix1-nstr1-nuni1-1))/(one-qs1)
        enddo
      case default
        call mpistop("unknown stretch type")
      end select
 
      select case(strtype2)
      case(0) ! uniform
        dxif2(:)=(xifmax2-xifmin2)/nxif2
        do ix2=1,nxif2
          xif2(ix2)=xifmin2+dxif2(ix2)*(ix2-half)
        enddo
      case(1) ! uni stretch
        qs2=qs2**(one/2**(refine_max_level-1))
        dxfirst2=(xifmax2-xifmin2)*(one-qs2)/(one-qs2**nxif2)
        dxif2(1)=dxfirst2
        do ix2=2,nxif1
          dxif2(ix2)=dxfirst2*qs2**(ix2-1)
          xif2(ix2)=dxif2(1)/(one-qs2)*(one-qs2**(ix2-1))+half*dxif2(ix2)
        enddo
      case(2) ! symm stretch
        ! base level, nbb = nstr + nuni + nstr
        nstr2=nstrb2*bnx2/2
        nuni2=nbb2-nstrb2*bnx2
        lenstr2=(xifmax2-xifmin2)/(2.d0+nuni2*(one-qs2)/(one-qs2**nstr2))
        dxfirst2=(xifmax2-xifmin2)/(dble(nuni2)+2.d0/(one-qs2)*(one-qs2**nstr2))
        dxmid2=dxfirst2
        ! refine_max level, numXI = nstr + nuni + nstr
        nstr2=nstr2*2**(refine_max_level-1)
        nuni2=nuni2*2**(refine_max_level-1)
        qs2=qs2**(one/2**(refine_max_level-1))
        dxfirst2=lenstr2*(one-qs2)/(one-qs2**nstr2)
        dxmid2=dxmid2/2**(refine_max_level-1)
        ! uniform center
        if(nuni2 .gt. 0) then
          do ix2=nstr2+1,nstr2+nuni2
            dxif2(ix2)=dxmid2
            xif2(ix2)=lenstr2+(dble(ix2)-0.5d0-nstr2)*dxif2(ix2)+xifmin2
          enddo
        endif
        ! left half
        do ix2=nstr2,1,-1
          dxif2(ix2)=dxfirst2*qs2**(nstr2-ix2)
          xif2(ix2)=xifmin2+lenstr2-dxif2(ix2)*half-dxfirst2*(one-qs2**(nstr2-ix2))/(one-qs2)
        enddo
        ! right half
        do ix2=nstr2+nuni2+1,nxif2
          dxif2(ix2)=dxfirst2*qs2**(ix2-nstr2-nuni2-1)
          xif2(ix2)=xifmax2-lenstr2+dxif2(ix2)*half+dxfirst2*(one-qs2**(ix2-nstr2-nuni2-1))/(one-qs2)
        enddo
      case default
        call mpistop("unknown stretch type")
      end select
 
      ! integrate EUV flux and get cell average flux for image
      if (datatype=='image_euv') then
        if (si_unit) then
          unitv=unit_velocity/1.0e3 ! km/s
        else
          unitv=unit_velocity/1.0e5 ! km/s
        endif
        numwi=2
        allocate(wi(nxif1,nxif2,numwi))
        allocate(euvs(nxif1,nxif2),euv(nxif1,nxif2))
        allocate(dpl(nxif1,nxif2),dpls(nxif1,nxif2))
        euvs=0.0d0
        euv=0.0d0
        dpl=0.d0
        dpls=0.d0
        do iigrid=1,igridstail; igrid=igrids(iigrid);
          call integrate_euv_datresol(igrid,nxif1,nxif2,xif1,xif2,dxif1,dxif2,fl,euvs,dpls)
        enddo
        numsi=nxif1*nxif2
        call mpi_allreduce(euvs,euv,numsi,mpi_double_precision, &
                           mpi_sum,icomm,ierrmpi)
        call mpi_allreduce(dpls,dpl,numsi,mpi_double_precision, &
                           mpi_sum,icomm,ierrmpi)
        do ix1=1,nxif1
          do ix2=1,nxif2
            if (euv(ix1,ix2)<smalldouble) euv(ix1,ix2)=zero
            if(euv(ix1,ix2)/=0) then
              dpl(ix1,ix2)=(dpl(ix1,ix2)/euv(ix1,ix2))*unitv
            else
              dpl(ix1,ix2)=0.d0
            endif
            if (abs(dpl(ix1,ix2))<smalldouble) dpl(ix1,ix2)=zero
          enddo
        enddo
        wi(:,:,1)=euv(:,:)
        wi(:,:,2)=dpl(:,:)
 
        call output_data(qunit,xif1,xif2,dxif1,dxif2,wi,nxif1,nxif2,numwi,datatype)
        deallocate(wi,euv,euvs,dpl,dpls)
      endif
 
      ! integrate SXR flux and get cell average flux for image
      if (datatype=='image_sxr') then
        if (si_unit) then
          arcsec=7.25d5
        else
          arcsec=7.25d7
        endif
        rhessi_rsl=2.3d0/instrument_resolution_factor
        numwi=1
        allocate(wi(nxif1,nxif2,numwi))
        allocate(sxrs(nxif1,nxif2),sxr(nxif1,nxif2))
        sxrs=0.0d0
        sxr=0.0d0
        do iigrid=1,igridstail; igrid=igrids(iigrid);
          call integrate_sxr_datresol(igrid,nxif1,nxif2,xif1,xif2,dxif1,dxif2,fl,sxrs)
        enddo
        numsi=nxif1*nxif2
        call mpi_allreduce(sxrs,sxr,numsi,mpi_double_precision, &
                           mpi_sum,icomm,ierrmpi)
 
        sxr=sxr*(rhessi_rsl*arcsec)**2 ! photons cm^-2 s^-1 pixel^-1
        do ix1=1,nxif1
          do ix2=1,nxif2
            if (sxr(ix1,ix2)<smalldouble) sxr(ix1,ix2)=zero
          enddo
        enddo
        wi(:,:,1)=sxr(:,:)
 
        call output_data(qunit,xif1,xif2,dxif1,dxif2,wi,nxif1,nxif2,numwi,datatype)
        deallocate(wi,sxr,sxrs)
      endif
 
      deallocate(xif1,xif2,dxif1,dxif2)
 
    end subroutine get_image_datresol
 
    subroutine integrate_sxr_datresol(igrid,nXIF1,nXIF2,xIF1,xIF2,dxIF1,dxIF2,fl,SXR)
      use mod_global_parameters
 
      integer, intent(in) :: igrid,nXIF1,nXIF2
      double precision, intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
      double precision, intent(in) :: dxIF1(nXIF1),dxIF2(nXIF2)
      type(te_fluid), intent(in) :: fl
      double precision, intent(out) :: SXR(nXIF1,nXIF2)
 
      integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
      double precision :: xb^L,xd^D
      double precision, allocatable :: flux(:^D&)
      double precision, allocatable :: dxb1(:^D&),dxb2(:^D&),dxb3(:^D&)
      double precision, allocatable :: SXRg(:,:),xg1(:),xg2(:),dxg1(:),dxg2(:)
      integer :: levelg,nXg1,nXg2,iXgmin1,iXgmax1,iXgmin2,iXgmax2,rft,iXg^D
      double precision :: SXRt,xc^L,xg^L,r2,area_1AU
      integer :: ixP^L,ixP^D
      integer :: direction_LOS
 
      if (los_phi==0 .and. los_theta==90) then
        direction_los=1
      else if (los_phi==90 .and. los_theta==90) then
        direction_los=2
      else
        direction_los=3
      endif
 
      ^d&ixomin^d=ixmlo^d\
      ^d&ixomax^d=ixmhi^d\
      ^d&iximin^d=ixglo^d\
      ^d&iximax^d=ixghi^d\
      ^d&xbmin^d=rnode(rpxmin^d_,igrid)\
      ^d&xbmax^d=rnode(rpxmax^d_,igrid)\
 
      allocate(flux(ixi^s))
      allocate(dxb1(ixi^s),dxb2(ixi^s),dxb3(ixi^s))
      dxb1(ixo^s)=ps(igrid)%dx(ixo^s,1)
      dxb2(ixo^s)=ps(igrid)%dx(ixo^s,2)
      dxb3(ixo^s)=ps(igrid)%dx(ixo^s,3)
      ! get local SXR flux
      call get_sxr(ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux,emin_sxr,emax_sxr)
 
      ! grid parameters
      levelg=ps(igrid)%level
      rft=2**(refine_max_level-levelg)
 
      ! fine table for storing EUV flux of current grid
      select case(direction_los)
      case(1)
        nxg1=iximax2*rft
        nxg2=iximax3*rft
      case(2)
        nxg1=iximax3*rft
        nxg2=iximax1*rft
      case(3)
        nxg1=iximax1*rft
        nxg2=iximax2*rft
      end select
      allocate(sxrg(nxg1,nxg2),xg1(nxg1),xg2(nxg2),dxg1(nxg1),dxg2(nxg2))
      sxrg=zero
      xg1=zero
      xg2=zero
 
      ! integrate for different direction
      select case(direction_los)
      case(1)
        do ix2=ixomin2,ixomax2
          ixgmin1=(ix2-1)*rft+1
          ixgmax1=ix2*rft
          do ix3=ixomin3,ixomax3
            ixgmin2=(ix3-1)*rft+1
            ixgmax2=ix3*rft
            sxrt=0.d0
            do ix1=ixomin1,ixomax1
              sxrt=sxrt+flux(ix^d)*dxb1(ix^d)*unit_length
            enddo
            sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=sxrt
          enddo
        enddo
      case(2)
        do ix3=ixomin3,ixomax3
          ixgmin1=(ix3-1)*rft+1
          ixgmax1=ix3*rft
          do ix1=ixomin1,ixomax1
            ixgmin2=(ix1-1)*rft+1
            ixgmax2=ix1*rft
            sxrt=0.d0
            do ix2=ixomin2,ixomax2
              sxrt=sxrt+flux(ix^d)*dxb2(ix^d)*unit_length
            enddo
            sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=sxrt
          enddo
        enddo
      case(3)
        do ix1=ixomin1,ixomax1
          ixgmin1=(ix1-1)*rft+1
          ixgmax1=ix1*rft
          do ix2=ixomin2,ixomax2
            ixgmin2=(ix2-1)*rft+1
            ixgmax2=ix2*rft
            sxrt=0.d0
            do ix3=ixomin3,ixomax3
              sxrt=sxrt+flux(ix^d)*dxb3(ix^d)*unit_length
            enddo
            sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=sxrt
          enddo
        enddo
      end select
 
      area_1au=2.81d27
      sxrg=sxrg/area_1au
 
      ! mapping grid data to global table 
      ! index ranges in local table
      select case(direction_los)
      case(1)
        ixgmin1=(ixomin2-1)*rft+1
        ixgmax1=ixomax2*rft
        ixgmin2=(ixomin3-1)*rft+1
        ixgmax2=ixomax3*rft
      case(2)
        ixgmin1=(ixomin3-1)*rft+1
        ixgmax1=ixomax3*rft
        ixgmin2=(ixomin1-1)*rft+1
        ixgmax2=ixomax1*rft
      case(3)
        ixgmin1=(ixomin1-1)*rft+1
        ixgmax1=ixomax1*rft
        ixgmin2=(ixomin2-1)*rft+1
        ixgmax2=ixomax2*rft
      end select
      ! index ranges in global table & mapping
      select case(direction_los)
      case(1)
        ixpmin1=(node(pig2_,igrid)-1)*rft*block_nx2+1
        ixpmax1=node(pig2_,igrid)*rft*block_nx2
        ixpmin2=(node(pig3_,igrid)-1)*rft*block_nx3+1
        ixpmax2=node(pig3_,igrid)*rft*block_nx3
      case(2)
        ixpmin1=(node(pig3_,igrid)-1)*rft*block_nx3+1
        ixpmax1=node(pig3_,igrid)*rft*block_nx3
        ixpmin2=(node(pig1_,igrid)-1)*rft*block_nx1+1
        ixpmax2=node(pig1_,igrid)*rft*block_nx1
      case(3)
        ixpmin1=(node(pig1_,igrid)-1)*rft*block_nx1+1
        ixpmax1=node(pig1_,igrid)*rft*block_nx1
        ixpmin2=(node(pig2_,igrid)-1)*rft*block_nx2+1
        ixpmax2=node(pig2_,igrid)*rft*block_nx2
      end select
      xg1(ixgmin1:ixgmax1)=xif1(ixpmin1:ixpmax1)
      xg2(ixgmin2:ixgmax2)=xif2(ixpmin2:ixpmax2)
      dxg1(ixgmin1:ixgmax1)=dxif1(ixpmin1:ixpmax1)
      dxg2(ixgmin2:ixgmax2)=dxif2(ixpmin2:ixpmax2)
      sxr(ixpmin1:ixpmax1,ixpmin2:ixpmax2)=sxr(ixpmin1:ixpmax1,ixpmin2:ixpmax2)+&
                                           sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)
 
      deallocate(flux,dxb1,dxb2,dxb3,sxrg,xg1,xg2,dxg1,dxg2)
 
    end subroutine integrate_sxr_datresol
 
    subroutine integrate_euv_datresol(igrid,nXIF1,nXIF2,xIF1,xIF2,dxIF1,dxIF2,fl,EUV,Dpl)
      use mod_global_parameters
 
      integer, intent(in) :: igrid,nXIF1,nXIF2
      double precision, intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
      double precision, intent(in) :: dxIF1(nXIF1),dxIF2(nXIF2)
      type(te_fluid), intent(in) :: fl
      double precision, intent(out) :: EUV(nXIF1,nXIF2),Dpl(nXIF1,nXIF2)
 
      integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
      double precision :: xb^L,xd^D
      double precision, allocatable :: flux(:^D&),v(:^D&),rho(:^D&)
      double precision, allocatable :: dxb1(:^D&),dxb2(:^D&),dxb3(:^D&)
      double precision, allocatable :: EUVg(:,:),Fvg(:,:),xg1(:),xg2(:),dxg1(:),dxg2(:)
      integer :: levelg,nXg1,nXg2,iXgmin1,iXgmax1,iXgmin2,iXgmax2,rft,iXg^D
      double precision :: EUVt,Fvt,xc^L,xg^L,r2
      integer :: ixP^L,ixP^D
      integer :: direction_LOS
 
      if (los_phi==0 .and. los_theta==90) then
        direction_los=1
      else if (los_phi==90 .and. los_theta==90) then
        direction_los=2
      else
        direction_los=3
      endif
 
      ^d&ixomin^d=ixmlo^d\
      ^d&ixomax^d=ixmhi^d\
      ^d&iximin^d=ixglo^d\
      ^d&iximax^d=ixghi^d\
      ^d&xbmin^d=rnode(rpxmin^d_,igrid)\
      ^d&xbmax^d=rnode(rpxmax^d_,igrid)\
 
      allocate(flux(ixi^s),v(ixi^s),rho(ixi^s))
      allocate(dxb1(ixi^s),dxb2(ixi^s),dxb3(ixi^s))
      dxb1(ixo^s)=ps(igrid)%dx(ixo^s,1)
      dxb2(ixo^s)=ps(igrid)%dx(ixo^s,2)
      dxb3(ixo^s)=ps(igrid)%dx(ixo^s,3)
      ! get local EUV flux and velocity
      call get_euv(wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
      flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2   ! adjust flux due to artifical change of resolution
      call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
      v(ixo^s)=-ps(igrid)%w(ixo^s,iw_mom(direction_los))/rho(ixo^s)
      deallocate(rho)
 
      ! grid parameters
      levelg=ps(igrid)%level
      rft=2**(refine_max_level-levelg)
 
      ! fine table for storing EUV flux of current grid
      select case(direction_los)
      case(1)
        nxg1=iximax2*rft
        nxg2=iximax3*rft
      case(2)
        nxg1=iximax3*rft
        nxg2=iximax1*rft
      case(3)
        nxg1=iximax1*rft
        nxg2=iximax2*rft
      end select
      allocate(euvg(nxg1,nxg2),fvg(nxg1,nxg2),xg1(nxg1),xg2(nxg2),dxg1(nxg1),dxg2(nxg2))
      euvg=zero
      fvg=zero
      xg1=zero
      xg2=zero
 
      ! integrate for different direction
      select case(direction_los)
      case(1)
        do ix2=ixomin2,ixomax2
          ixgmin1=(ix2-1)*rft+1
          ixgmax1=ix2*rft
          do ix3=ixomin3,ixomax3
            ixgmin2=(ix3-1)*rft+1
            ixgmax2=ix3*rft
            euvt=0.d0
            fvt=0.d0
            do ix1=ixomin1,ixomax1
              euvt=euvt+flux(ix^d)*dxb1(ix^d)*unit_length
              fvt=fvt+flux(ix^d)*dxb1(ix^d)*unit_length*v(ix^d)
            enddo
            euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=euvt
            fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=fvt
          enddo
        enddo
      case(2)
        do ix3=ixomin3,ixomax3
          ixgmin1=(ix3-1)*rft+1
          ixgmax1=ix3*rft
          do ix1=ixomin1,ixomax1
            ixgmin2=(ix1-1)*rft+1
            ixgmax2=ix1*rft
            euvt=0.d0
            fvt=0.d0
            do ix2=ixomin2,ixomax2
              euvt=euvt+flux(ix^d)*dxb2(ix^d)*unit_length
              fvt=fvt+flux(ix^d)*dxb2(ix^d)*unit_length*v(ix^d)
            enddo
            euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=euvt
            fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=fvt
          enddo
        enddo
      case(3)
        do ix1=ixomin1,ixomax1
          ixgmin1=(ix1-1)*rft+1
          ixgmax1=ix1*rft
          do ix2=ixomin2,ixomax2
            ixgmin2=(ix2-1)*rft+1
            ixgmax2=ix2*rft
            euvt=0.d0
            fvt=0.d0
            do ix3=ixomin3,ixomax3
              euvt=euvt+flux(ix^d)*dxb3(ix^d)*unit_length
              fvt=fvt+flux(ix^d)*dxb3(ix^d)*unit_length*v(ix^d)
            enddo
            euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=euvt
            fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=fvt
          enddo
        enddo
      end select
      if (si_unit) then
        euvg=euvg*1.d2
        fvg=fvg*1.d2
      endif
 
      ! mapping grid data to global table 
      ! index ranges in local table
      select case(direction_los)
      case(1)
        ixgmin1=(ixomin2-1)*rft+1
        ixgmax1=ixomax2*rft
        ixgmin2=(ixomin3-1)*rft+1
        ixgmax2=ixomax3*rft
      case(2)
        ixgmin1=(ixomin3-1)*rft+1
        ixgmax1=ixomax3*rft
        ixgmin2=(ixomin1-1)*rft+1
        ixgmax2=ixomax1*rft
      case(3)
        ixgmin1=(ixomin1-1)*rft+1
        ixgmax1=ixomax1*rft
        ixgmin2=(ixomin2-1)*rft+1
        ixgmax2=ixomax2*rft
      end select
      ! index ranges in global table & mapping
      select case(direction_los)
      case(1)
        ixpmin1=(node(pig2_,igrid)-1)*rft*block_nx2+1
        ixpmax1=node(pig2_,igrid)*rft*block_nx2
        ixpmin2=(node(pig3_,igrid)-1)*rft*block_nx3+1
        ixpmax2=node(pig3_,igrid)*rft*block_nx3
      case(2)
        ixpmin1=(node(pig3_,igrid)-1)*rft*block_nx3+1
        ixpmax1=node(pig3_,igrid)*rft*block_nx3
        ixpmin2=(node(pig1_,igrid)-1)*rft*block_nx1+1
        ixpmax2=node(pig1_,igrid)*rft*block_nx1
      case(3)
        ixpmin1=(node(pig1_,igrid)-1)*rft*block_nx1+1
        ixpmax1=node(pig1_,igrid)*rft*block_nx1
        ixpmin2=(node(pig2_,igrid)-1)*rft*block_nx2+1
        ixpmax2=node(pig2_,igrid)*rft*block_nx2
      end select
      xg1(ixgmin1:ixgmax1)=xif1(ixpmin1:ixpmax1)
      xg2(ixgmin2:ixgmax2)=xif2(ixpmin2:ixpmax2)
      dxg1(ixgmin1:ixgmax1)=dxif1(ixpmin1:ixpmax1)
      dxg2(ixgmin2:ixgmax2)=dxif2(ixpmin2:ixpmax2)
      euv(ixpmin1:ixpmax1,ixpmin2:ixpmax2)=euv(ixpmin1:ixpmax1,ixpmin2:ixpmax2)+&
                                           euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)
      dpl(ixpmin1:ixpmax1,ixpmin2:ixpmax2)=dpl(ixpmin1:ixpmax1,ixpmin2:ixpmax2)+&
                                           fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)
 
      deallocate(flux,v,dxb1,dxb2,dxb3,euvg,fvg,xg1,xg2,dxg1,dxg2)
 
    end subroutine integrate_euv_datresol
 
    subroutine get_image(qunit,datatype,fl)
      ! integrate emission flux along line of sight (LOS) 
      ! in a 3D simulation box and get a 2D EUV image
      use mod_global_parameters
      use mod_constants
 
      integer, intent(in) :: qunit
      type(te_fluid), intent(in) :: fl
      character(20), intent(in) :: datatype
 
      integer :: ix^D,numXI1,numXI2,numWI
      double precision :: xImin1,xImax1,xImin2,xImax2,xIcent1,xIcent2,dxI
      double precision, allocatable :: xI1(:),xI2(:),dxI1(:),dxI2(:)
      double precision, allocatable :: wI(:,:,:),wIs(:,:,:),EM(:,:),WLB(:,:,:)
      double precision :: vec_temp1(1:3),vec_temp2(1:3)
      double precision :: vec_z(1:3),vec_cor(1:3),xI_cor(1:2)
      double precision :: res,LOS_psi,r_max,r_loc
 
      integer :: mass
      character (30) :: ion
      double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
      double precision :: arcsec,RHESSI_rsl,LASCO_rsl,pixel,R_occult,smallflux
      integer :: iigrid,igrid,i,j,numSI,iw
      logical :: emit
 
      if (coordinate==spherical) then
        call init_vectors_spherical()
      else
        ! cartesian
        call init_vectors_cartesian()
      endif
 
      ! calculate domain of the image
      if (coordinate==spherical) then
        ximin1=-abs(xprobmax1)
        ximin2=-abs(xprobmax1)
        ximax1=abs(xprobmax1)
        ximax2=abs(xprobmax1)
      else
        ! calculate domain of the image
        do ix1=1,2
          if (ix1==1) vec_cor(1)=xprobmin1
          if (ix1==2) vec_cor(1)=xprobmax1
          do ix2=1,2
            if (ix2==1) vec_cor(2)=xprobmin2
            if (ix2==2) vec_cor(2)=xprobmax2
            do ix3=1,2
              if (ix3==1) vec_cor(3)=xprobmin3
              if (ix3==2) vec_cor(3)=xprobmax3
              if (big_image) then
                r_loc=(vec_cor(1)-x_origin(1))**2
                r_loc=r_loc+(vec_cor(2)-x_origin(2))**2
                r_loc=r_loc+(vec_cor(3)-x_origin(3))**2
                r_loc=sqrt(r_loc)
                if (ix1==1 .and. ix2==1 .and. ix3==1) then
                  r_max=r_loc
                else
                  r_max=max(r_max,r_loc)
                endif
              else
                call get_cor_image(vec_cor,xi_cor)
                if (ix1==1 .and. ix2==1 .and. ix3==1) then
                  ximin1=xi_cor(1)
                  ximax1=xi_cor(1)
                  ximin2=xi_cor(2)
                  ximax2=xi_cor(2)
                else
                  ximin1=min(ximin1,xi_cor(1))
                  ximax1=max(ximax1,xi_cor(1))
                  ximin2=min(ximin2,xi_cor(2))
                  ximax2=max(ximax2,xi_cor(2))
                endif
              endif
            enddo
          enddo
        enddo
        if (big_image) then
          ximin1=-r_max
          ximin2=-r_max
          ximax1=r_max
          ximax2=r_max
        endif
      endif
      xicent1=(ximin1+ximax1)/2.d0
      xicent2=(ximin2+ximax2)/2.d0
 
      ! tables for image
      if (si_unit) then
        arcsec=7.25d5/unit_length
      else
        arcsec=7.25d7/unit_length
      endif
      if (datatype=='image_euv') then
        call get_line_info(wavelength,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
        dxi=spacersl*arcsec  ! intrument resolution of image
        smallflux=smalldouble
      else if (datatype=='image_sxr') then
        rhessi_rsl=2.3d0/instrument_resolution_factor
        dxi=rhessi_rsl*arcsec
        smallflux=1.d-40
      else if (datatype=='image_whitelight') then
        if (whitelight_instrument=='LASCO/C1') then
          lasco_rsl=5.6d0/instrument_resolution_factor
          r_occult=1.1d0
        else if (whitelight_instrument=='LASCO/C2') then
          lasco_rsl=11.4d0/instrument_resolution_factor
          r_occult=2.d0
        else if (whitelight_instrument=='LASCO/C3') then
          lasco_rsl=56.d0/instrument_resolution_factor
          r_occult=3.7d0
        else
          call mpistop('Whitelight synthesis: instrument is not supported!')
        endif
        dxi=lasco_rsl*arcsec
        if (r_occultor>1.d0) r_occult=r_occultor
        r_occult=r_occult*const_rsun/unit_length
        smallflux=1.d-20
      endif
      numxi1=8*ceiling((ximax1-xicent1)/dxi/8.d0)
      ximin1=xicent1-numxi1*dxi
      ximax1=xicent1+numxi1*dxi
      numxi1=numxi1*2
      numxi2=8*ceiling((ximax2-xicent2)/dxi/8.d0)
      ximin2=xicent2-numxi2*dxi
      ximax2=xicent2+numxi2*dxi
      numxi2=numxi2*2
      allocate(xi1(numxi1),xi2(numxi2),dxi1(numxi1),dxi2(numxi2))
      do ix1=1,numxi1
        xi1(ix1)=ximin1+dxi*(ix1-half)
        dxi1(ix1)=dxi
      enddo
      do ix2=1,numxi2
        xi2(ix2)=ximin2+dxi*(ix2-half)
        dxi2(ix2)=dxi
      enddo
 
      ! calculate emission
      if (datatype=='image_euv' .or. datatype=='image_sxr') then
        numwi=1
        allocate(wi(numxi1,numxi2,numwi),wis(numxi1,numxi2,numwi),em(numxi1,numxi2))
        wi=zero
        wis=zero
        em=zero
        if (coordinate==cartesian) then
          do iigrid=1,igridstail; igrid=igrids(iigrid);
            call integrate_emission_cartesian(igrid,numxi1,numxi2,xi1,xi2,dxi,fl,datatype,em)
          enddo
        else
          do iigrid=1,igridstail; igrid=igrids(iigrid);
            call integrate_emission_spherical(igrid,numxi1,numxi2,xi1,xi2,dxi,fl,datatype,em)
          enddo
        endif
        do ix1=1,numxi1
          do ix2=1,numxi2
            if (em(ix1,ix2)>smallflux) wis(ix1,ix2,1)=em(ix1,ix2)
          enddo
        enddo
        numsi=numxi1*numxi2*numwi
        call mpi_allreduce(wis,wi,numsi,mpi_double_precision,mpi_sum,icomm,ierrmpi)
        if (activate_unit_arcsec) then
          xi1=xi1/arcsec
          dxi1=dxi1/arcsec
          xi2=xi2/arcsec
          dxi2=dxi2/arcsec
        endif
        call output_data(qunit,xi1,xi2,dxi1,dxi2,wi,numxi1,numxi2,numwi,datatype)
        deallocate(wi,wis,em)
      else if (datatype=='image_whitelight') then
        numwi=2
        allocate(wi(numxi1,numxi2,numwi),wis(numxi1,numxi2,numwi),wlb(numxi1,numxi2,numwi))
        wi=zero
        wis=zero
        wlb=zero
        if (coordinate==spherical) then
          do iigrid=1,igridstail; igrid=igrids(iigrid);
            call integrate_whitelight_spherical(igrid,numxi1,numxi2,numwi,xi1,xi2,dxi,fl,datatype,wlb)
          enddo
        endif
        do ix1=1,numxi1
          do ix2=1,numxi2
            if (wlb(ix1,ix2,1)>smallflux) then 
              wis(ix1,ix2,1)=wlb(ix1,ix2,1)
              wis(ix1,ix2,2)=wlb(ix1,ix2,2)
            endif
          enddo
        enddo
        numsi=numxi1*numxi2*numwi
        call mpi_allreduce(wis,wi,numsi,mpi_double_precision,mpi_sum,icomm,ierrmpi)
        if (activate_unit_arcsec) then
          xi1=xi1/arcsec
          dxi1=dxi1/arcsec
          xi2=xi2/arcsec
          dxi2=dxi2/arcsec
        endif
        call output_data(qunit,xi1,xi2,dxi1,dxi2,wi,numxi1,numxi2,numwi,datatype)
        deallocate(wi,wis,wlb)
      endif
 
      deallocate(xi1,xi2,dxi1,dxi2)
 
    end subroutine get_image
 
    subroutine integrate_emission_cartesian(igrid,numXI1,numXI2,xI1,xI2,dxI,fl,datatype,EM)
      integer, intent(in) :: igrid,numXI1,numXI2
      double precision, intent(in) :: xI1(numXI1),xI2(numXI2)
      double precision, intent(in) :: dxI
      type(te_fluid), intent(in) :: fl
      character(20), intent(in) :: datatype
      double precision, intent(inout) :: EM(numXI1,numXI2)
 
      integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
      double precision :: xb^L,xd^D
      double precision, allocatable :: flux(:^D&)
      double precision :: res
      integer :: ixP^L,ixP^D,nSubC^D,iSubC^D
      double precision :: xSubP1,xSubP2,dxSubP,xerf^L,fluxsubC
      double precision :: xSubC(1:3),dxSubC^D,xCent(1:2)
 
      integer :: mass
      double precision :: logTe
      character (30) :: ion
      double precision :: lineCent
      double precision :: sigma_PSF,spaceRsl,wlRsl,sigma0,factor,wslit
      double precision :: arcsec,pixel,RHESSI_rsl,area_1AU
      double precision :: aa,bb
 
      ^d&ixomin^d=ixmlo^d\
      ^d&ixomax^d=ixmhi^d\
      ^d&iximin^d=ixglo^d\
      ^d&iximax^d=ixghi^d\
      ^d&xbmin^d=rnode(rpxmin^d_,igrid)\
      ^d&xbmax^d=rnode(rpxmax^d_,igrid)\
 
      if (si_unit) then
        arcsec=7.25d5/unit_length
      else
        arcsec=7.25d7/unit_length
      endif
 
      allocate(flux(ixi^s))
      if (datatype=='image_euv') then
        ! get local EUV flux and velocity
        call get_euv(wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
        flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2   ! adjust flux due to artifical change of resolution
        call get_line_info(wavelength,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
        pixel=spacersl*arcsec
        sigma0=sigma_psf*pixel
      else if (datatype=='image_sxr') then
        ! get local SXR flux photons cm^-3 s^-1 (cgs) or photons m^-3 s^-1 (SI)
        call get_sxr(ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux,emin_sxr,emax_sxr) 
        rhessi_rsl=2.3d0/instrument_resolution_factor
        sigma_psf=1.d0
        pixel=rhessi_rsl*arcsec
        sigma0=sigma_psf*pixel
        area_1au=2.81d27
      endif
 
      ! integrate emission
      {do ix^d=ixomin^d,ixomax^d\}
        ^d&nsubc^d=1;
        ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(vec_xi1(^d))/(dxi/2.d0)));
        ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(vec_xi2(^d))/(dxi/2.d0)));
        ^d&dxsubc^d=ps(igrid)%dx(ix^dd,^d)/nsubc^d;
        if (datatype=='image_euv') then
          if (si_unit) then
            fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length*1.d2/dxi/dxi  ! DN s^-1
          else
            fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length/dxi/dxi  ! DN s^-1
          endif
        else if (datatype=='image_sxr') then
          ! sub-cell SXR flux at 1 AU [photons s^-1 cm^-2]
          fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length**3/area_1au
        endif
        if (fluxsubc>smalldouble) then
          ! dividing a cell to several parts to get more accurate integrating values
          {do isubc^d=1,nsubc^d\}
            ^d&xsubc(^d)=ps(igrid)%x(ix^dd,^d)-half*ps(igrid)%dx(ix^dd,^d)+(isubc^d-half)*dxsubc^d;
            ! mapping the 3D coordinate to location at the image
            call get_cor_image(xsubc,xcent)
            ! distribution at nearby pixels
            ixp1=floor((xcent(1)-(xi1(1)-half*dxi))/dxi)+1
            ixp2=floor((xcent(2)-(xi2(1)-half*dxi))/dxi)+1
            ixpmin1=max(1,ixp1-3)
            ixpmax1=min(ixp1+3,numxi1)
            ixpmin2=max(1,ixp2-3)
            ixpmax2=min(ixp2+3,numxi2)
            do ixp1=ixpmin1,ixpmax1
              do ixp2=ixpmin2,ixpmax2
                xerfmin1=((xi1(ixp1)-half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0) 
                xerfmax1=((xi1(ixp1)+half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
                xerfmin2=((xi2(ixp2)-half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
                xerfmax2=((xi2(ixp2)+half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
                factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
                em(ixp1,ixp2)=em(ixp1,ixp2)+fluxsubc*factor
              enddo !ixP2
            enddo !ixP1
          {enddo\} !iSubC
        endif
      {enddo\} !ix
 
      deallocate(flux)
    end subroutine integrate_emission_cartesian
 
    subroutine integrate_emission_spherical(igrid,numXI1,numXI2,xI1,xI2,dxI,fl,datatype,EM)
      integer, intent(in) :: igrid,numXI1,numXI2
      double precision, intent(in) :: xI1(numXI1),xI2(numXI2)
      double precision, intent(in) :: dxI
      type(te_fluid), intent(in) :: fl
      character(20), intent(in) :: datatype
      double precision, intent(inout) :: EM(numXI1,numXI2)
 
      integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
      double precision, allocatable :: flux(:^D&),Ne(:^D&)
      integer :: ixP^L,ixP^D,nSubC^D,iSubC^D
      double precision :: xSubP1,xSubP2,dxSubP,xerf^L,fluxsubC,RsubC
      double precision :: TBsubC,PBsubC
      double precision :: xSubC(1:3),dxSubC^D,xCent(1:2),xSubC_car(1:3)
      double precision :: R_thick,dotp,dvolume,R_occult,Rc
      double precision :: dxl(1:3),x_sph(1:3),dx_sph(1:3)
      double precision :: unitv_r(1:3),unitv_theta(1:3),unitv_phi(1:3)
      logical :: sun_back_side,emit
 
      integer :: mass
      double precision :: logTe
      character (30) :: ion
      double precision :: lineCent
      double precision :: sigma_PSF,spaceRsl,wlRsl,sigma0,factor,wslit
      double precision :: RHESSI_rsl,area_1AU,arcsec,pixel
 
      ^d&ixomin^d=ixmlo^d;
      ^d&ixomax^d=ixmhi^d;
      ^d&iximin^d=ixglo^d;
      ^d&iximax^d=ixghi^d;
 
      if (si_unit) then
        arcsec=7.25d5/unit_length
      else
        arcsec=7.25d7/unit_length
      endif
 
      allocate(flux(ixi^s))
      if (datatype=='image_euv') then
        ! get local EUV flux and velocity
        call get_euv(wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
        flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2   ! adjust flux due to artifical change of resolution
        call get_line_info(wavelength,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
        pixel=spacersl*arcsec
        sigma0=sigma_psf*pixel
      else if (datatype=='image_sxr') then
        ! get local SXR flux photons cm^-3 s^-1 (cgs) or photons m^-3 s^-1 (SI)
        call get_sxr(ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux,emin_sxr,emax_sxr) 
        rhessi_rsl=2.3d0/instrument_resolution_factor
        sigma_psf=1.d0
        pixel=rhessi_rsl*arcsec
        sigma0=sigma_psf*pixel
        area_1au=2.81d27
      endif
 
      ! integrate emission
      r_thick=r_opt_thick*const_rsun/unit_length
      {do ix^d=ixomin^d,ixomax^d\}
        x_sph(1:3)=ps(igrid)%x(ix^d,1:3)
        dx_sph(1:3)=ps(igrid)%dx(ix^d,1:3)
        dxl(1)=dx_sph(1) ! cell size in length
        dxl(2)=x_sph(1)*dx_sph(2) ! cell size in length
        dxl(3)=x_sph(1)*dsin(x_sph(2))*dx_sph(3) ! cell size in length
        ! dividing a cell to several sub-cells to get more accurate integrating values
        ^d&nsubc^d=1;
        call get_unit_vector_spherical(x_sph,unitv_r,unitv_theta,unitv_phi)
        call dot_product_loc(unitv_r,vec_xi1,dotp)
        nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
        call dot_product_loc(unitv_r,vec_xi2,dotp)
        nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
        call dot_product_loc(unitv_theta,vec_xi1,dotp)
        nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
        call dot_product_loc(unitv_theta,vec_xi2,dotp)
        nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
        call dot_product_loc(unitv_phi,vec_xi1,dotp)
        nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
        call dot_product_loc(unitv_phi,vec_xi2,dotp)
        nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
 
        ! integrate sub-cells
        do isubc1=1,nsubc1
          ! sub-cell center coordinate in spherical
          xsubc(1)=x_sph(1)-half*dx_sph(1)+(isubc1-half)*dx_sph(1)/nsubc1
          rsubc=xsubc(1)
          dxsubc1=dx_sph(1)/nsubc1 ! sub-cell size in length
          do isubc2=1,nsubc2
            ! sub-cell center coordinate in spherical
            xsubc(2)=x_sph(2)-half*dx_sph(2)+(isubc2-half)*dx_sph(2)/nsubc2
            dxsubc2=xsubc(1)*dx_sph(2)/nsubc2 ! sub-cell size in length
            dxsubc3=xsubc(1)*dsin(xsubc(2))*dx_sph(3)/nsubc3  ! sub-cell size in length
            dvolume=dxsubc1*dxsubc2*dxsubc3
            if (datatype=='image_euv') then
              if (si_unit) then
                fluxsubc=flux(ix^d)*dvolume*unit_length*1.d2/dxi/dxi  ! DN s^-1
              else
                fluxsubc=flux(ix^d)*dvolume*unit_length/dxi/dxi  ! DN s^-1
              endif
            else if (datatype=='image_sxr') then
              ! sub-cell SXR flux at 1 AU [photons s^-1 cm^-2]
              fluxsubc=flux(ix^d)*dvolume*unit_length**3/area_1au
            endif
            ! enter integration if flux large enough
            if (fluxsubc>smalldouble) then
              do isubc3=1,nsubc3
                ! sub-cell center coordinate in spherical
                xsubc(3)=x_sph(3)-half*dx_sph(3)+(isubc3-half)*dx_sph(3)/nsubc3
                call get_cor_image_spherical(xsubc,xcent)
                rc=dsqrt(xcent(1)**2+xcent(2)**2) ! distance to sun center (on the image plane)
                !
                ! whether the local emitted photons can arrive the telescope
                call spherical_to_cartesian(xsubc,xsubc_car)
                call dot_product_loc(vec_los,xsubc_car,dotp)
                sun_back_side=.true.
                if (dotp<0.d0) sun_back_side=.false.
                ! whether the local emission can reach the telescope
                if (sun_back_side) then
                  emit=.false.
                  if (rc>r_thick) emit=.true.
                else
                  emit=.true.
                  if (xsubc(1)<=r_thick) emit=.false.
                endif
                !
                if (emit) then
                  ! mapping the 3D coordinate to location at the image
                  ! distribution at nearby pixels
                  ixp1=floor((xcent(1)-(xi1(1)-half*dxi))/dxi)+1
                  ixp2=floor((xcent(2)-(xi2(1)-half*dxi))/dxi)+1
                  ixpmin1=max(1,ixp1-3)
                  ixpmax1=min(ixp1+3,numxi1)
                  ixpmin2=max(1,ixp2-3)
                  ixpmax2=min(ixp2+3,numxi2)
                  do ixp1=ixpmin1,ixpmax1
                    do ixp2=ixpmin2,ixpmax2
                      xerfmin1=((xi1(ixp1)-half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
                      xerfmax1=((xi1(ixp1)+half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
                      xerfmin2=((xi2(ixp2)-half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
                      xerfmax2=((xi2(ixp2)+half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
                      factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
                      em(ixp1,ixp2)=em(ixp1,ixp2)+fluxsubc*factor
                    enddo !ixP2
                  enddo !ixP1
                endif !emit
              enddo !iSubC3
            endif !smallflux
          enddo !iSubC2
        enddo !iSubC1
      {enddo\} !ix
 
      deallocate(flux)
 
    end subroutine integrate_emission_spherical
 
    subroutine integrate_whitelight_spherical(igrid,numXI1,numXI2,numWI,xI1,xI2,dxI,fl,datatype,WLB)
      integer, intent(in) :: igrid,numXI1,numXI2,numWI
      double precision, intent(in) :: xI1(numXI1),xI2(numXI2)
      double precision, intent(in) :: dxI
      type(te_fluid), intent(in) :: fl
      character(20), intent(in) :: datatype
      double precision, intent(inout) :: WLB(numXI1,numXI2,numWI)
 
      integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
      double precision, allocatable :: flux(:^D&),Ne(:^D&)
      integer :: ixP^L,ixP^D,nSubC^D,iSubC^D
      double precision :: xSubP1,xSubP2,dxSubP,xerf^L,fluxsubC,RsubC
      double precision :: sigma_PSF,sigma0,arcsec,pixel,LASCO_rsl
      double precision :: A,B,C,D,Rc,Ne0,TBsubC,PBsubC,factor
      double precision :: R_thick,dotp,dvolume,R_occult
      double precision :: xSubC(1:3),dxSubC^D,xCent(1:2),xSubC_car(1:3)
      double precision :: dxl(1:3),x_sph(1:3),dx_sph(1:3)
      double precision :: unitv_r(1:3),unitv_theta(1:3),unitv_phi(1:3)
      logical :: emit
 
      ^d&ixomin^d=ixmlo^d;
      ^d&ixomax^d=ixmhi^d;
      ^d&iximin^d=ixglo^d;
      ^d&iximax^d=ixghi^d;
 
      if (si_unit) then
        arcsec=7.25d5/unit_length
      else
        arcsec=7.25d7/unit_length
      endif
 
      allocate(ne(ixi^s))
      if (whitelight_instrument=='LASCO/C1') then
        lasco_rsl=5.6d0/instrument_resolution_factor
        r_occult=1.1d0
      else if (whitelight_instrument=='LASCO/C2') then
        lasco_rsl=11.4d0/instrument_resolution_factor
        r_occult=2.d0
      else if (whitelight_instrument=='LASCO/C3') then
        lasco_rsl=56.d0/instrument_resolution_factor
        r_occult=3.7d0
      endif
      if (r_occultor>1.d0) r_occult=r_occultor
      r_occult=r_occult*const_rsun/unit_length
      call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,ne)
      sigma_psf=1.d0
      pixel=lasco_rsl*arcsec
      sigma0=sigma_psf*pixel
 
      ! integrate emission
      r_thick=r_opt_thick*const_rsun/unit_length
      {do ix^d=ixomin^d,ixomax^d\}
        x_sph(1:3)=ps(igrid)%x(ix^d,1:3)
        dx_sph(1:3)=ps(igrid)%dx(ix^d,1:3)
        dxl(1)=dx_sph(1) ! cell size in length
        dxl(2)=x_sph(1)*dx_sph(2) ! cell size in length
        dxl(3)=x_sph(1)*dsin(x_sph(2))*dx_sph(3) ! cell size in length
        ne0=ne(ix^d)*unit_numberdensity
        ! dividing a cell to several sub-cells to get more accurate integrating values
        ^d&nsubc^d=1;
        call get_unit_vector_spherical(x_sph,unitv_r,unitv_theta,unitv_phi)
        call dot_product_loc(unitv_r,vec_xi1,dotp)
        nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
        call dot_product_loc(unitv_r,vec_xi2,dotp)
        nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
        call dot_product_loc(unitv_theta,vec_xi1,dotp)
        nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
        call dot_product_loc(unitv_theta,vec_xi2,dotp)
        nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
        call dot_product_loc(unitv_phi,vec_xi1,dotp)
        nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
        call dot_product_loc(unitv_phi,vec_xi2,dotp)
        nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
 
        ! integrate sub-cells
        do isubc1=1,nsubc1
          ! sub-cell center coordinate in spherical
          xsubc(1)=x_sph(1)-half*dx_sph(1)+(isubc1-half)*dx_sph(1)/nsubc1
          rsubc=xsubc(1)
          dxsubc1=dx_sph(1)/nsubc1 ! sub-cell size in length
          call get_thomson_parameters(rsubc,a,b,c,d)
          do isubc2=1,nsubc2
            ! sub-cell center coordinate in spherical
            xsubc(2)=x_sph(2)-half*dx_sph(2)+(isubc2-half)*dx_sph(2)/nsubc2
            dxsubc2=xsubc(1)*dx_sph(2)/nsubc2 ! sub-cell size in length
            dxsubc3=xsubc(1)*dsin(xsubc(2))*dx_sph(3)/nsubc3  ! sub-cell size in length
            dvolume=dxsubc1*dxsubc2*dxsubc3
            do isubc3=1,nsubc3
              ! sub-cell center coordinate in spherical
              xsubc(3)=x_sph(3)-half*dx_sph(3)+(isubc3-half)*dx_sph(3)/nsubc3
              call get_cor_image_spherical(xsubc,xcent)
              rc=dsqrt(xcent(1)**2+xcent(2)**2) ! distance to sun center (on the image plane)
              ! whether the local emitted photons can arrive the telescope
              emit=.false.
              if (rc>r_occult) then 
                emit=.true.
                ! scaterring flux from cm^-3 of plasma
                call get_whitelight_thomson(rsubc,rc,ne0,a,b,c,d,tbsubc,pbsubc)
                tbsubc=tbsubc*dvolume*unit_length/dxi/dxi
                pbsubc=pbsubc*dvolume*unit_length/dxi/dxi
                if (tbsubc<1.d-20) emit=.false.
              endif
              if (emit) then
                ! mapping the 3D coordinate to location at the image
                ! distribution at nearby pixels
                ixp1=floor((xcent(1)-(xi1(1)-half*dxi))/dxi)+1
                ixp2=floor((xcent(2)-(xi2(1)-half*dxi))/dxi)+1
                ixpmin1=max(1,ixp1-3)
                ixpmax1=min(ixp1+3,numxi1)
                ixpmin2=max(1,ixp2-3)
                ixpmax2=min(ixp2+3,numxi2)
                do ixp1=ixpmin1,ixpmax1
                  do ixp2=ixpmin2,ixpmax2
                    xerfmin1=((xi1(ixp1)-half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
                    xerfmax1=((xi1(ixp1)+half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
                    xerfmin2=((xi2(ixp2)-half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
                    xerfmax2=((xi2(ixp2)+half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
                    factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
                    wlb(ixp1,ixp2,1)=wlb(ixp1,ixp2,1)+tbsubc*factor
                    wlb(ixp1,ixp2,2)=wlb(ixp1,ixp2,2)+pbsubc*factor
                  enddo !ixP2
                enddo !ixP1
              endif
            enddo !iSubC3
          enddo !iSubC2
        enddo !iSubC1
      {enddo\} !ix
 
      deallocate(ne)
 
    end subroutine integrate_whitelight_spherical
 
    subroutine get_thomson_parameters(Rl,A,B,C,D)
      ! parameters given in Billings 1968
      use mod_constants
      double precision, intent(in) :: Rl
      double precision, intent(inout) :: A,B,C,D
 
      double precision :: sinO,cosO,sinO2,cosO2,tmp
 
      sino=const_rsun/(rl*unit_length)
      sino2=sino**2
      coso2=1.d0-sino2
      coso=abs(dsqrt(coso2))
      tmp=log((1.d0+sino)/coso)
      a=coso*sino2
      b=-(1.d0-3.d0*sino2-(coso2/sino)*(1.d0+3.d0*sino2)*tmp)/8.d0
      c=4.d0/3.d0-coso-coso*coso2/3.d0
      d=(5.d0+sino2-(coso2/sino)*(5.d0-sino2)*tmp)/8.d0
 
    end subroutine get_thomson_parameters
 
    subroutine get_whitelight_thomson(Rl,Rin,Ne,A,B,C,D,fluxTB,fluxPB)
      ! use the method in SSW/eltheory
      double precision, intent(in) :: Rl,Rin,Ne,A,B,C,D
      double precision, intent(inout) :: fluxTB,fluxPB
 
      double precision :: const,u,Bt,Br,PB,TB,sinchi2
 
      u=0.63d0
      const=1.24878d-25/(1.d0-u/3.d0)
      sinchi2=(rin/rl)**2
      bt=const*(c+u*(d-c))
      pb=const*sinchi2*((a+u*(b-a)))
      br=bt-pb
      tb=bt+br
      fluxtb=tb*ne
      fluxpb=pb*ne
 
    end subroutine get_whitelight_thomson
 
    subroutine get_unit_vector_spherical(x_sph,unitv_r,unitv_theta,unitv_phi)
      double precision, intent(in) :: x_sph(1:3)
      double precision, intent(inout) :: unitv_r(1:3),unitv_theta(1:3),unitv_phi(1:3)
 
      unitv_r(1)=dsin(x_sph(2))*dcos(x_sph(3))
      unitv_r(2)=dsin(x_sph(2))*dsin(x_sph(3))
      unitv_r(3)=dcos(x_sph(2))
      unitv_theta(1)=dcos(x_sph(2))*dcos(x_sph(3))
      unitv_theta(2)=dcos(x_sph(2))*dsin(x_sph(3))
      unitv_theta(3)=-dsin(x_sph(2))
      unitv_phi(1)=-dsin(x_sph(3))
      unitv_phi(2)=dcos(x_sph(3))
      unitv_phi(3)=zero
 
    end subroutine get_unit_vector_spherical
 
    subroutine output_data(qunit,xO1,xO2,dxO1,dxO2,wO,nXO1,nXO2,nWO,datatype)
      ! change the format of data and write data
      use mod_global_parameters
 
      integer, intent(in) :: qunit,nXO1,nXO2,nWO
      double precision, intent(in) :: dxO1(nxO1),dxO2(nxO2)
      double precision, intent(in) :: xO1(nXO1),xO2(nxO2)
      double precision, intent(inout) :: wO(nXO1,nXO2,nWO)
      character(20), intent(in) :: datatype
 
      integer :: nPiece,nP1,nP2,nC1,nC2,nWC
      integer :: piece_nmax1,piece_nmax2,ix1,ix2,j,ipc,ixc1,ixc2
      double precision, allocatable :: xC(:,:,:,:),wC(:,:,:,:),dxC(:,:,:,:)
 
      ! clean small values
      do ix1=1,nxo1
        do ix2=1,nxo2
          do j=1,nwo
            if (abs(wo(ix1,ix2,j))<smalldouble) wo(ix1,ix2,j)=zero
          enddo
        enddo
      enddo
 
      ! how many cells in each grid
      if (dat_resolution) then
        if (datatype=='image_euv' .or. datatype=='image_sxr') then
          if (los_phi==0 .and. los_theta==90) then
            piece_nmax1=block_nx2
            piece_nmax2=block_nx3
          else if (los_phi==90 .and. los_theta==90) then
            piece_nmax1=block_nx3
            piece_nmax2=block_nx1
          else
            piece_nmax1=block_nx1
            piece_nmax2=block_nx2
          endif
        else if (datatype=='spectrum_euv') then
          piece_nmax1=16
          if (direction_slit==1) then
            piece_nmax2=block_nx1
          else if (direction_slit==2) then
            piece_nmax2=block_nx2
          else
            piece_nmax2=block_nx3
          endif
        endif
      else
        piece_nmax1=20
        piece_nmax2=20
      endif
      loopn1: do j=piece_nmax1,1,-1
        if(mod(nxo1,j)==0) then
          nc1=j
          exit loopn1
        endif
      enddo loopn1
      loopn2: do j=piece_nmax2,1,-1
        if(mod(nxo2,j)==0) then
          nc2=j
          exit loopn2
        endif
      enddo loopn2
      ! how many grids
      np1=nxo1/nc1
      np2=nxo2/nc2
      npiece=np1*np2
      nwc=nwo
 
      ! output images
      select case(convert_type)
        case('EIvtuCCmpi','ESvtuCCmpi','SIvtuCCmpi','WIvtuCCmpi')
          ! put data into grids
          allocate(xc(npiece,nc1,nc2,2))
          allocate(dxc(npiece,nc1,nc2,2))
          allocate(wc(npiece,nc1,nc2,nwo))
          do ipc=1,npiece
            do ixc1=1,nc1
              do ixc2=1,nc2
                ix1=mod(ipc-1,np1)*nc1+ixc1
                ix2=floor(1.0*(ipc-1)/np1)*nc2+ixc2
                xc(ipc,ixc1,ixc2,1)=xo1(ix1)
                xc(ipc,ixc1,ixc2,2)=xo2(ix2)
                dxc(ipc,ixc1,ixc2,1)=dxo1(ix1)
                dxc(ipc,ixc1,ixc2,2)=dxo2(ix2)
                do j=1,nwc
                  wc(ipc,ixc1,ixc2,j)=wo(ix1,ix2,j)
                enddo
              enddo
            enddo
          enddo
          ! write data into vtu file
          call write_image_vtucc(qunit,xc,wc,dxc,npiece,nc1,nc2,nwc,datatype)
          deallocate(xc,dxc,wc)
        case('EIvtiCCmpi','ESvtiCCmpi','SIvtiCCmpi','WIvtiCCmpi')
          if (sum(stretch_type(:))>0 .and. dat_resolution) then
            call mpistop("Error in synthesize emission: vti is not supported for dat resolution")
          else
            call write_image_vticc(qunit,xo1,xo2,dxo1,dxo2,wo,nxo1,nxo2,nwo,nc1,nc2)
          endif
        case default
          call mpistop("Error in synthesize emission: Unknown convert_type")
      end select
 
    end subroutine output_data
  }
 
    subroutine write_image_vticc(qunit,xO1,xO2,dxO1,dxO2,wO,nXO1,nXO2,nWO,nC1,nC2)
      ! write image data to vti
      use mod_global_parameters
 
      integer, intent(in) :: qunit,nXO1,nXO2,nWO,nC1,nC2
      double precision, intent(in) :: xO1(nXO1),xO2(nxO2)
      double precision, intent(in) :: dxO1(nxO1),dxO2(nxO2)
      double precision, intent(in) :: wO(nXO1,nXO2,nWO)
 
      double precision :: origin(1:3), spacing(1:3)
      integer :: wholeExtent(1:6),extent(1:6)
      integer :: nP1,nP2,iP1,iP2,iw
      integer :: ixC1,ixC2,ixCmin1,ixCmax1,ixCmin2,ixCmax2
 
      integer :: filenr
      logical :: fileopen
      character (70) :: subname,wname,vname,nameL,nameS
      character (len=std_len) :: filename
      integer :: mass
      double precision :: logTe
      character (30) :: ion
      double precision :: line_center
      double precision :: spatial_rsl,spectral_rsl,sigma_PSF,wslit
 
 
      origin(1)=xo1(1)-0.5d0*dxo1(1)
      origin(2)=xo2(1)-0.5d0*dxo2(1)
      origin(3)=zero
      spacing(1)=dxo1(1)
      spacing(2)=dxo2(1)
      spacing(3)=zero
      wholeextent=0
      wholeextent(2)=nxo1
      wholeextent(4)=nxo2
      np1=nxo1/nc1
      np2=nxo2/nc2
 
      ! get information of emission line
      if (convert_type=='EIvtiCCmpi') then
        call get_line_info(wavelength,ion,mass,logte,line_center,spatial_rsl,spectral_rsl,sigma_psf,wslit)
      else if (convert_type=='ESvtiCCmpi') then
        call get_line_info(spectrum_wl,ion,mass,logte,line_center,spatial_rsl,spectral_rsl,sigma_psf,wslit)
      endif
 
      if (mype==0) then
        inquire(qunit,opened=fileopen)
        if(.not.fileopen)then
          ! generate filename 
          filenr=snapshotini
          if (autoconvert) filenr=snapshotnext
          if (convert_type=='EIvtiCCmpi') then
            write(filename,'(a,i4.4,a)') trim(filename_euv),filenr,".vti"
          else if (convert_type=='SIvtiCCmpi') then
            write(filename,'(a,i4.4,a)') trim(filename_sxr),filenr,".vti"
          else if (convert_type=='WIvtiCCmpi') then
            write(filename,'(a,i4.4,a)') trim(filename_whitelight),filenr,".vti"
          else if (convert_type=='ESvtiCCmpi') then
            write(filename,'(a,i4.4,a)') trim(filename_spectrum),filenr,".vti"
          endif
          open(qunit,file=filename,status='unknown',form='formatted')
        endif
 
        ! generate xml header
        write(qunit,'(a)')'<?xml version="1.0"?>'
        write(qunit,'(a)',advance='no') '<VTKFile type="ImageData"'
        write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
        write(qunit,'(a,3(1pe14.6),a,6(i10),a,3(1pe14.6),a)')'  <ImageData Origin="',&
              origin,'" WholeExtent="',wholeextent,'" Spacing="',spacing,'">'
        ! file info        
        write(qunit,'(a)')'<FieldData>'
        write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
                           'NumberOfTuples="1" format="ascii">'
        write(qunit,*) real(global_time*time_convert_factor)
        write(qunit,'(a)')'</DataArray>'
        if (convert_type=='EIvtiCCmpi' .or. convert_type=='ESvtiCCmpi') then
          write(qunit,'(2a)')'<DataArray type="Float32" Name="logT" ',&
                             'NumberOfTuples="1" format="ascii">'
          write(qunit,*) real(logte)
          write(qunit,'(a)')'</DataArray>'
        endif
        write(qunit,'(a)')'</FieldData>'
        ! pixel/cell data
        do ip1=1,np1
          do ip2=1,np2
            extent=0
            extent(1)=(ip1-1)*nc1
            extent(2)=ip1*nc1
            extent(3)=(ip2-1)*nc2
            extent(4)=ip2*nc2
            ixcmin1=extent(1)+1
            ixcmax1=extent(2)
            ixcmin2=extent(3)+1
            ixcmax2=extent(4)
            write(qunit,'(a,6(i10),a)') &
                  '<Piece Extent="',extent,'">'
            write(qunit,'(a)')'<CellData>'
            do iw=1,nwo
              ! variable name
              if (convert_type=='EIvtiCCmpi') then
                if (wavelength<100) then
                  write(vname,'(a,i2)') "AIA",wavelength
                else if (wavelength<1000) then
                  write(vname,'(a,i3)') "AIA",wavelength
                else
                  write(vname,'(a,i4)') "IRIS",wavelength
                endif
                if (iw==2 .and. dat_resolution) vname='Doppler_velocity'
              else if (convert_type=='SIvtiCCmpi') then
                if (emin_sxr<10 .and. emax_sxr<10) then
                  write(vname,'(a,i1,a,i1,a)') "SXR",emin_sxr,"-",emax_sxr,"keV"
                else if (emin_sxr<10 .and. emax_sxr>=10) then
                  write(vname,'(a,i1,a,i2,a)') "SXR",emin_sxr,"-",emax_sxr,"keV"
                else
                  write(vname,'(a,i2,a,i2,a)') "SXR",emin_sxr,"-",emax_sxr,"keV"
                endif
              else if (convert_type=='WIvtiCCmpi') then
                if (iw==1) write(vname,'(a)')'B'
                if (iw==2) write(vname,'(a)')'pB'
              else if (convert_type=='ESvtiCCmpi') then
                if (spectrum_wl==1354) then
                  write(vname,'(a,i4)') "SG",spectrum_wl
                else
                  write(vname,'(a,i3)') "EIS",spectrum_wl
                endif
              endif
              write(qunit,'(a,a,a)')&
                '<DataArray type="Float64" Name="',trim(vname),'" format="ascii">'
              write(qunit,'(200(1pe14.6))') ((wo(ixc1,ixc2,iw),ixc1=ixcmin1,ixcmax1),ixc2=ixcmin2,ixcmax2)
              write(qunit,'(a)')'</DataArray>'
            enddo
            write(qunit,'(a)')'</CellData>'
            write(qunit,'(a)')'</Piece>'
          enddo
        enddo
        ! end
        write(qunit,'(a)')'</ImageData>'
        write(qunit,'(a)')'</VTKFile>'
        close(qunit)
      endif
 
    end subroutine write_image_vticc
 
    subroutine write_image_vtucc(qunit,xC,wC,dxC,nPiece,nC1,nC2,nWC,datatype)
      ! write image data to vtu
      use mod_global_parameters
 
      integer, intent(in) :: qunit
      integer, intent(in) :: nPiece,nC1,nC2,nWC
      double precision, intent(in) :: xC(nPiece,nC1,nC2,2),dxC(nPiece,nc1,nc2,2)
      double precision, intent(in) :: wC(nPiece,nC1,nC2,nWC)
      character(20), intent(in) :: datatype
 
      integer :: nP1,nP2
      double precision :: xP(nPiece,nC1+1,nC2+1,2)
      integer :: filenr
      logical :: fileopen
      character (70) :: subname,wname,vname,nameL,nameS
      character (len=std_len) :: filename
      integer :: ixC1,ixC2,ixP,ix1,ix2,j
      integer :: nc,np,icel,VTK_type
      integer :: mass
      double precision :: logTe
      character (30) :: ion
      double precision :: line_center
      double precision :: spatial_rsl,spectral_rsl,sigma_PSF,wslit
 
      np1=nc1+1
      np2=nc2+1
      np=np1*np2
      nc=nc1*nc2
      ! cell corner location     
      do ixp=1,npiece
        do ix1=1,np1
          do ix2=1,np2
            if (ix1<np1) xp(ixp,ix1,ix2,1)=xc(ixp,ix1,1,1)-0.5d0*dxc(ixp,ix1,1,1)
            if (ix1==np1) xp(ixp,ix1,ix2,1)=xc(ixp,ix1-1,1,1)+0.5d0*dxc(ixp,ix1-1,1,1)
            if (ix2<np2) xp(ixp,ix1,ix2,2)=xc(ixp,1,ix2,2)-0.5d0*dxc(ixp,1,ix2,2)
            if (ix2==np2) xp(ixp,ix1,ix2,2)=xc(ixp,1,ix2-1,2)+0.5d0*dxc(ixp,1,ix2-1,2)
          enddo
        enddo
      enddo
      ! get information of emission line
      if (datatype=='image_euv') then
        call get_line_info(wavelength,ion,mass,logte,line_center,spatial_rsl,spectral_rsl,sigma_psf,wslit)
      else if (datatype=='spectrum_euv') then
        call get_line_info(spectrum_wl,ion,mass,logte,line_center,spatial_rsl,spectral_rsl,sigma_psf,wslit)
      endif
 
      if (mype==0) then
        inquire(qunit,opened=fileopen)
        if(.not.fileopen)then
          ! generate filename 
          filenr=snapshotini
          if (autoconvert) filenr=snapshotnext
          if (datatype=='image_euv') then
            write(filename,'(a,i4.4,a)') trim(filename_euv),filenr,".vtu"
          else if (datatype=='image_sxr') then
            write(filename,'(a,i4.4,a)') trim(filename_sxr),filenr,".vtu"
          else if (datatype=='image_whitelight') then
            write(filename,'(a,i4.4,a)') trim(filename_whitelight),filenr,".vtu"
          else if (datatype=='spectrum_euv') then
            write(filename,'(a,i4.4,a)') trim(filename_spectrum),filenr,".vtu"
          endif
          open(qunit,file=filename,status='unknown',form='formatted')
        endif
        ! generate xml header
        write(qunit,'(a)')'<?xml version="1.0"?>'
        write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
        write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
        write(qunit,'(a)')'<UnstructuredGrid>'
        write(qunit,'(a)')'<FieldData>'
        write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
                           'NumberOfTuples="1" format="ascii">'
        write(qunit,*) real(global_time*time_convert_factor)
        write(qunit,'(a)')'</DataArray>'
        if (datatype=='image_euv' .or. datatype=='spectrum_euv') then
          write(qunit,'(2a)')'<DataArray type="Float32" Name="logT" ',&
                             'NumberOfTuples="1" format="ascii">'
          write(qunit,*) real(logte)
          write(qunit,'(a)')'</DataArray>'
        endif
        write(qunit,'(a)')'</FieldData>'
        do ixp=1,npiece
          write(qunit,'(a,i7,a,i7,a)') &
                '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
          write(qunit,'(a)')'<CellData>'
          do j=1,nwc
            if (datatype=='image_euv') then
              if (j==1) then
                if (wavelength<100) then
                  write(vname,'(a,i2)') "AIA",wavelength
                else if (wavelength<1000) then
                  write(vname,'(a,i3)') "AIA",wavelength
                else
                  write(vname,'(a,i4)') "IRIS",wavelength
                endif
              endif
              if (j==2 .and. dat_resolution) vname='Doppler_velocity'
            else if (datatype=='image_sxr') then
              if (emin_sxr<10 .and. emax_sxr<10) then
                write(vname,'(a,i1,a,i1,a)') "SXR",emin_sxr,"-",emax_sxr,"keV"
              else if (emin_sxr<10 .and. emax_sxr>=10) then
                write(vname,'(a,i1,a,i2,a)') "SXR",emin_sxr,"-",emax_sxr,"keV"
              else
                write(vname,'(a,i2,a,i2,a)') "SXR",emin_sxr,"-",emax_sxr,"keV"
              endif
            else if (datatype=='image_whitelight') then
              write(vname,'(a)')'whitelight'
            else if (datatype=='spectrum_euv') then
              if (spectrum_wl==1354) then
                write(vname,'(a,i4)') "SG",spectrum_wl
              else
                write(vname,'(a,i3)') "EIS",spectrum_wl
              endif
            endif
            write(qunit,'(a,a,a)')&
              '<DataArray type="Float64" Name="',trim(vname),'" format="ascii">'
            write(qunit,'(200(1pe14.6))') ((wc(ixp,ixc1,ixc2,j),ixc1=1,nc1),ixc2=1,nc2)
            write(qunit,'(a)')'</DataArray>'
          enddo
          write(qunit,'(a)')'</CellData>'
          write(qunit,'(a)')'<Points>'
          write(qunit,'(a)')'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
          do ix2=1,np2
            do ix1=1,np1
              if (datatype=='image_euv' .and. dat_resolution) then
                if (los_phi==0 .and. los_theta==90) then
                  write(qunit,'(3(1pe14.6))') 0.d0,xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2)
                else if (los_phi==90 .and. los_theta==90) then
                  write(qunit,'(3(1pe14.6))') xp(ixp,ix1,ix2,2),0.d0,xp(ixp,ix1,ix2,1)
                else
                  write(qunit,'(3(1pe14.6))') xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2),0.d0
                endif
              else if (datatype=='image_sxr' .and. dat_resolution) then
                if (los_phi==0 .and. los_theta==90) then
                  write(qunit,'(3(1pe14.6))') 0.d0,xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2)
                else if (los_phi==90 .and. los_theta==90) then
                  write(qunit,'(3(1pe14.6))') xp(ixp,ix1,ix2,2),0.d0,xp(ixp,ix1,ix2,1)
                else
                  write(qunit,'(3(1pe14.6))') xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2),0.d0
                endif
              else
                write(qunit,'(3(1pe14.6))') xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2),0.d0
              endif
            enddo
          enddo
          write(qunit,'(a)')'</DataArray>'
          write(qunit,'(a)')'</Points>'
          ! connetivity part
          write(qunit,'(a)')'<Cells>'
          write(qunit,'(a)')'<DataArray type="Int32" Name="connectivity" format="ascii">'
          do ix2=1,nc2
            do ix1=1,nc1
              write(qunit,'(4(i7))') ix1-1+(ix2-1)*np1,ix1+(ix2-1)*np1,&
                                     ix1-1+ix2*np1,ix1+ix2*np1
            enddo
          enddo
          write(qunit,'(a)')'</DataArray>'
          ! offsets data array
          write(qunit,'(a)')'<DataArray type="Int32" Name="offsets" format="ascii">'
          do icel=1,nc
            write(qunit,'(i7)') icel*(2**2)
          enddo
          write(qunit,'(a)')'</DataArray>'
          ! VTK cell type data array
          write(qunit,'(a)')'<DataArray type="Int32" Name="types" format="ascii">'
          ! VTK_LINE=3; VTK_PIXEL=8; VTK_VOXEL=11 -> vtk-syntax
          vtk_type=8
          do icel=1,nc
            write(qunit,'(i2)') vtk_type
          enddo
          write(qunit,'(a)')'</DataArray>'
          write(qunit,'(a)')'</Cells>'
          write(qunit,'(a)')'</Piece>'
        enddo
        write(qunit,'(a)')'</UnstructuredGrid>'
        write(qunit,'(a)')'</VTKFile>'
        close(qunit)
      endif
    end subroutine write_image_vtucc
 
    subroutine dot_product_loc(vec1,vec2,res)
      double precision, intent(in) :: vec1(1:3),vec2(1:3)
      double precision, intent(out) :: res
 
      res=vec1(1)*vec2(1)+vec1(2)*vec2(2)+vec1(3)*vec2(3)
 
    end subroutine dot_product_loc
 
    subroutine cross_product_loc(vec_in1,vec_in2,vec_out)
      double precision, intent(in) :: vec_in1(1:3),vec_in2(1:3)
      double precision, intent(out) :: vec_out(1:3)
 
      vec_out(1)=vec_in1(2)*vec_in2(3)-vec_in1(3)*vec_in2(2)
      vec_out(2)=vec_in1(3)*vec_in2(1)-vec_in1(1)*vec_in2(3)
      vec_out(3)=vec_in1(1)*vec_in2(2)-vec_in1(2)*vec_in2(1)
 
    end subroutine cross_product_loc
 
    subroutine init_vectors_spherical()
      integer :: j
      double precision :: LOS_psi
      double precision :: vec_car(1:3),vec_z(1:3),vec_temp1(1:3),vec_temp2(1:3)
      double precision :: vec_LOS_sph(1:3),vec_xI1_sph(1:3),vec_xI2_sph(1:3)
 
      ! antiparallel to LOS in spherical
      vec_los(1)=1.d0
      vec_los(2)=dpi*los_theta/180.d0
      vec_los(3)=dpi*los_phi/180.d0
      ! LOS in cartesian
      call spherical_to_cartesian(vec_los,vec_car)
      vec_los=-vec_car
 
      ! theta=0 in cartesian
      vec_z(:)=zero
      vec_z(3)=1.d0
 
      ! x direction for image
      if (los_theta==zero) then
        vec_temp1(1)=1.d0
        vec_temp1(2)=dpi/2.d0
        vec_temp1(3)=dpi*los_phi/180.d0
        call spherical_to_cartesian(vec_temp1,vec_car)
        vec_temp1=-vec_car
        call cross_product_loc(vec_temp1,vec_z,vec_xi1)
      else
        call cross_product_loc(vec_los,vec_z,vec_xi1)
      endif
 
      ! y direction for image
      call cross_product_loc(vec_xi1,vec_los,vec_xi2)
 
      ! rotate the image
      vec_temp1=vec_xi1/sqrt(vec_xi1(1)**2+vec_xi1(2)**2+vec_xi1(3)**2)
      vec_temp2=vec_xi2/sqrt(vec_xi2(1)**2+vec_xi2(2)**2+vec_xi2(3)**2)
      los_psi=dpi*image_rotate/180.d0
      vec_xi1=vec_temp1*cos(los_psi)-vec_temp2*sin(los_psi)
      vec_xi2=vec_temp2*cos(los_psi)+vec_temp1*sin(los_psi)
 
      do j=1,3
        if (abs(vec_los(j))<smalldouble) vec_los(j)=zero
        if (abs(vec_xi1(j))<smalldouble) vec_xi1(j)=zero
        if (abs(vec_xi2(j))<smalldouble) vec_xi2(j)=zero
      enddo
 
      call cartesian_to_spherical(vec_los,vec_los_sph)
      call cartesian_to_spherical(vec_xi1,vec_xi1_sph)
      call cartesian_to_spherical(vec_xi2,vec_xi2_sph)
      vec_los_sph(2:3)=vec_los_sph(2:3)*180.d0/dpi
      vec_xi1_sph(2:3)=vec_xi1_sph(2:3)*180.d0/dpi
      vec_xi2_sph(2:3)=vec_xi2_sph(2:3)*180.d0/dpi
 
      if (mype==0) write(*,'(a,f3.1,f6.1,f6.1,a)') ' LOS vector (spherial): [',vec_los_sph(1),vec_los_sph(2),vec_los_sph(3),']'
      if (mype==0) write(*,'(a,f3.1,f6.1,f6.1,a)') ' xI1 vector (spherial): [',vec_xi1_sph(1),vec_xi1_sph(2),vec_xi1_sph(3),']'
      if (mype==0) write(*,'(a,f3.1,f6.1,f6.1,a)') ' xI2 vector (spherial): [',vec_xi2_sph(1),vec_xi2_sph(2),vec_xi2_sph(3),']'
 
    end subroutine init_vectors_spherical
 
    subroutine spherical_to_cartesian(vec_sph,vec_car)
      ! angles in rad
      double precision, intent(in) :: vec_sph(1:3)
      double precision, intent(inout) :: vec_car(1:3)
 
      vec_car(1)=vec_sph(1)*dsin(vec_sph(2))*dcos(vec_sph(3))
      vec_car(2)=vec_sph(1)*dsin(vec_sph(2))*dsin(vec_sph(3))
      vec_car(3)=vec_sph(1)*dcos(vec_sph(2))
 
    end subroutine spherical_to_cartesian
 
    subroutine cartesian_to_spherical(vec_car,vec_sph)
      ! angles in rad
      double precision, intent(in) :: vec_car(1:3)
      double precision, intent(inout) :: vec_sph(1:3)
 
      vec_sph(1)=dsqrt(vec_car(1)**2+vec_car(2)**2+vec_car(3)**2)
      vec_sph(2)=dacos(vec_car(3)/vec_sph(1))
      vec_sph(3)=atan2(vec_car(2),vec_car(3))
 
    end subroutine cartesian_to_spherical
 
    subroutine init_vectors_cartesian()
      integer :: j
      double precision :: LOS_psi
      double precision :: vec_z(1:3),vec_temp1(1:3),vec_temp2(1:3)
 
      ! vectors for image coordinate
      vec_los(1)=-cos(dpi*los_phi/180.d0)*sin(dpi*los_theta/180.d0)
      vec_los(2)=-sin(dpi*los_phi/180.d0)*sin(dpi*los_theta/180.d0)
      vec_los(3)=-cos(dpi*los_theta/180.d0)
      do j=1,3
        if (abs(vec_los(j))<=smalldouble) vec_los(j)=zero
      enddo
      vec_z(:)=zero
      vec_z(3)=1.d0
      if (los_theta==zero) then
        vec_xi1(1)=cos(dpi*los_phi/180.d0)
        vec_xi1(2)=sin(dpi*los_phi/180.d0)
        vec_xi1(3)=zero
      else
        call cross_product_loc(vec_los,vec_z,vec_xi1)
      endif
      call cross_product_loc(vec_xi1,vec_los,vec_xi2)
      vec_temp1=vec_xi1/sqrt(vec_xi1(1)**2+vec_xi1(2)**2+vec_xi1(3)**2)
      vec_temp2=vec_xi2/sqrt(vec_xi2(1)**2+vec_xi2(2)**2+vec_xi2(3)**2)
      los_psi=dpi*image_rotate/180.d0
      vec_xi1=vec_temp1*cos(los_psi)-vec_temp2*sin(los_psi)
      vec_xi2=vec_temp2*cos(los_psi)+vec_temp1*sin(los_psi)
 
      do j=1,3
        if (abs(vec_xi1(j))<smalldouble) vec_xi1(j)=zero
        if (abs(vec_xi2(j))<smalldouble) vec_xi2(j)=zero
      enddo
 
      if (mype==0) write(*,'(a,f5.2,f6.2,f6.2,a)') ' LOS vector: [',vec_los(1),vec_los(2),vec_los(3),']'
      if (mype==0) write(*,'(a,f5.2,f6.2,f6.2,a)') ' xI1 vector: [',vec_xi1(1),vec_xi1(2),vec_xi1(3),']'
      if (mype==0) write(*,'(a,f5.2,f6.2,f6.2,a)') ' xI2 vector: [',vec_xi2(1),vec_xi2(2),vec_xi2(3),']'
 
    end subroutine init_vectors_cartesian
 
    subroutine get_cor_image_spherical(x_3D_sph,x_image)
      double precision, intent(in) :: x_3D_sph(1:3)
      double precision, intent(inout) :: x_image(1:2)
      double precision :: res,res_origin
      double precision :: x_3D(1:3)
 
      call spherical_to_cartesian(x_3d_sph,x_3d)
      call dot_product_loc(x_3d,vec_xi1,res)
      x_image(1)=res
      call dot_product_loc(x_3d,vec_xi2,res)
      x_image(2)=res
 
    end subroutine get_cor_image_spherical
 
    subroutine get_cor_image(x_3D,x_image)
      double precision, intent(in) :: x_3D(1:3)
      double precision, intent(inout) :: x_image(1:2)
      double precision :: res,res_origin
 
      call dot_product_loc(x_3d,vec_xi1,res)
      call dot_product_loc(x_origin,vec_xi1,res_origin)
      x_image(1)=res-res_origin
      call dot_product_loc(x_3d,vec_xi2,res)
      call dot_product_loc(x_origin,vec_xi2,res_origin)
      x_image(2)=res-res_origin
 
    end subroutine get_cor_image
 
end module mod_thermal_emission