source: Sophya/trunk/Cosmo/RadioBeam/plpkn.pic@ 3973

########################################################################
### Script de trace de PNoise(k) et reponse dans le plan (u,v) 
###  de l'interferometre a partir du fichier PPF produit par pknoise.cc
###    Janvier 2011 ,  BAORadio/Reza
########################################################################

echo ' -----> plpknew.pic : opening PPF files ... "
delobjs *
# openppf ../../PkNoise/cmvhpkz.ppf 
# openppf ../pkz0p25.ppf
openppf ../pkz0p7.ppf
rename hpkz hpkz0p7
openppf ../pkz1p0.ppf
rename hpkz hpkz1p0

defscript pknopen 
  set f $1
  set nm $2
  echo '------ pknopen File= ' $f ' name=' $nm 
  openppf $f
  rename h1dnoise h1dn$nm
  rename h1drep h1dr$nm
  rename fracmodok fracok$nm
  mv dtnoise pkn$nm
  rename h2drep h2dr$nm
endscript 

pknopen pknf11x11.ppf A121d 
pknopen pknnan128.ppf B128d 
pknopen pknconfC.ppf  C129d 
pknopen pknD50m.ppf   D50m 
pknopen pknD75m.ppf   D75m 
pknopen pknD100m.ppf   D100m 
pknopen pknD200m.ppf   D200m 
pknopen pknf20x20.ppf E400d
pknopen pknf4cyl.ppf F4cyl
pknopen pknf4cylp.ppf F4cylp
pknopen pknf8cyl.ppf G8cyl
pknopen pknf8cylp.ppf G8cylp


listobjs 

#  A z = 0.7
z = 0.7
c = 3.e5
H = 102
Da = 2488 
nu21 = 1.42e9

echo '----> Executing anapkn.pic'
exec anapkn.pic
setup5 
# scalewz 0.7 2488 102 1
# scalewz 0.25 989.13 80.26 1
scalewz 1 3300 120.5 1
Da = 3300 


y1 = 400*$cct21
y2 = 7e4*$cct21
xyl = "xylimits=0.005,0.5,$y1,$y2 logx logy minorticks"
# H = 71.9 km/s/Mpc = 1.0271 x 70 km/s
h70 = 1.0271
h70cube = $h70*$h70*$h70
set kk pow(10.,x)/$h70

defscript plpklss
  scalewz 1 3300 120.5 1
#  scalewz 0.7 2488 102 1
  convpk2t21
  cct21h70 = $cct21*$h70cube
  n/plot hpkz1p0.val*${cct21h70}%$kk ! ! "notit nsta connectpoints black $xyl line=solid,2 "
#  n/plot hpkz0p7.val*${cct21h70}%$kk ! ! "notit nsta connectpoints black $xyl line=solid,2 "

  addtext 0.01 4000 '[ P(k)-LSS z=1.0 ]' 'font=helvetica,bolditalic,16 black'
  setaxelabels 'k_comov  (h70 Mpc^-1)'  'P21(k) mK^2 x (Mpc/h70)^3' 'font=helvetica,bolditalic,16'

endscript  

setaxesatt 'minorticks font=helvetica,bold,16 autofontsize'

Rad2Deg = 180/3.141596
R2D2 = $Rad2Deg*$Rad2Deg
Tsys = 50 
Delnu = 1.e6

nu = $nu21/(1+$z)
pi23 = 8.*Pi*Pi*Pi

Lambda = 0.21*(1+$z)
Lam2 = $Lambda*$Lambda 

Dx = 100
Dy = 100

defscript calcul 
  Dxol = $Dx/$Lambda 
  Dyol = $Dy/$Lambda 

#   FOV = (1.2*1.2*$Lam2/$Dx/$Dy) 
  FOV = ($Lam2/$Dx/$Dy) 

  FOVDEG = $FOV*$R2D2 
  NPointing = 10000/$FOVDEG
  tinteg = 365*86400/$NPointing
  PNOISE = 2.*$Tsys*$Tsys/$tinteg/$Dxol/$Dyol
  PNOISE = $PNOISE*$Da*$Da*$c/$H*(1+$z)/$nu
#   PNOISE = $PNOISE*$Da*$Da*$c/$H*(1+$z)/$nu/$pi23
  PNOISE = $PNOISE*1.e6
  PNOISE = $PNOISE*1.05
  echo " FOV = $FOV deg^2  NPointing= $NPointing"
  echo " tinteg= $tinteg sec  PNOISE= $PNOISE mK^2"
endscript 


defscript plnoiseDishes
  Dx = 75
  Dy = 75
  calcul 
  plot2d pknD75m  k/$Da pnoise*$PNOISE/5 fracmodok>0.2 'same line=solid,2 cpts red nsta notit'
  Dx = 100
  Dy = 100
  calcul 
  plot2d pknD100m  k/$Da pnoise*$PNOISE/20 fracmodok>0.2 'same line=solid,2 cpts blueviolet nsta notit'

  plot2d pknD100m  k/$Da pnoise*$PNOISE/100 fracmodok>0.2 'same line=solid,2 cpts blue nsta notit'



  set lines ( '(a) 75m Dish, 5 beams'  '(b) 100m Dish, 10 beams'  '(b) 100m Dish, 100 beams'  )
  set cols ( red blueviolet blue )
  textdrawer lines cols 'frame font=helvetica,bold,16  inset=0.1,0.3,0.15,0.35' 
  settitle ' PNoise(k) : 5/10/100 beams/polar @z=1' ' ' 'font=helvetica,bold,16' 
endscript

defscript plnoiseEFGH
  Dx = 12*0.9
  Dy = 0.844248*0.8
  calcul 
  plot2d pknF4cyl k/$Da pnoise*$PNOISE fracmodok>0.2 'same line=solid,2 cpts forestgreen nsta notit' 
  plot2d pknF4cylp k/$Da pnoise*$PNOISE fracmodok>0.2 'same line=solid,2 cpts green nsta notit' 

  plot2d pknG8cyl k/$Da pnoise*$PNOISE fracmodok>0.2 'same line=solid,2 cpts violetred nsta notit' 
  plot2d pknG8cylp k/$Da pnoise*$PNOISE fracmodok>0.2 'same line=solid,2 cpts violet nsta notit' 

  Dx = 5*0.9
  Dy = 5*0.9
  calcul 
  plot2d pknE400d k/$Da pnoise*$PNOISE fracmodok>0.2 'same line=solid,2 cpts turquoise nsta notit'

  set lines (  '(e) 20x20:400xD=5m' '(f) 4Cyl-12mx85m, 400 rec/pol' '(fp) 4Cylp-12mx85m, 400 rec/pol'  '(g) 8Cyl-12mx105m, 960 rec/pol' '(gp) 8Cylp-12mx105m, 960 rec/pol' )
  set cols ( turquoise forestgreen green violetred violet  )
  textdrawer lines cols 'frame font=helvetica,bold,16  inset=0.1,0.3,0.15,0.35' 
  settitle ' PNoise(k) : Dishes/Cylinders, 400/400/960 recv/pol @z=1' ' ' 'font=helvetica,bold,16' 
endscript 

defscript plnoiseABCD
  Dx = 5*0.9
  Dy = 5*0.9
  calcul 
  plot2d pknA121d k/$Da pnoise*$PNOISE fracmodok>0.2 'same line=solid,2 cpts magenta nsta notit'
  plot2d pknB128d k/$Da pnoise*$PNOISE fracmodok>0.2 'same line=solid,2 cpts red nsta notit'
  plot2d pknC129d k/$Da pnoise*$PNOISE fracmodok>0.2 'same line=solid,2 cpts orange nsta notit'

  Dx = 75
  Dy = 75
  calcul 
  plot2d pknD75m  k/$Da pnoise*$PNOISE/100 fracmodok>0.2 'same line=solid,2 cpts blueviolet nsta notit'

  set lines (  '(a) 11x11:121xD=5m' '(b) 128xD=5m' '(c) 129xD=5m' '(d) 75m Dish, 100 beams' )
  set cols ( magenta red orange blueviolet )
  textdrawer lines cols 'frame font=helvetica,bold,16  inset=0.1,0.3,0.15,0.35' 
  settitle ' PNoise(k) : 121,128,129,100 receiver/polar @z=1' ' ' 'font=helvetica,bold,16' 
endscript

defscript plfracABCD
  xylf = "xylimits=0.005,0.5,0.,1. minorticks"
  plot2d pknA121d k/$Da  fracmodok fracmodok>0. '$xylf line=solid,2 cpts magenta nsta notit'
  plot2d pknB128d k/$Da  fracmodok fracmodok>0. 'same line=solid,2 cpts red nsta notit'
  plot2d pknC129d k/$Da  fracmodok fracmodok>0. 'same line=solid,2 cpts orange nsta notit'
  plot2d pknD75m  k/$Da  fracmodok fracmodok>0. 'same line=solid,2 cpts violetblue nsta notit'
endscript

defscript plfracEFGH
  xylf = "xylimits=0.005,0.5,0.,1. minorticks"
  plot2d pknF4cyl k/$Da  fracmodok fracmodok>0. '$xylf line=solid,2 cpts forestgreen nsta notit'
  plot2d pknF4cylp k/$Da  fracmodok fracmodok>0. 'same line=solid,2 cpts green nsta notit'
  plot2d pknG8cyl k/$Da  fracmodok fracmodok>0. 'same line=solid,2 cpts gold nsta notit'
  plot2d pknG8cylp  k/$Da  fracmodok fracmodok>0. 'same line=solid,2 cpts yellow nsta notit'
  plot2d pknE400d  k/$Da  fracmodok fracmodok>0. 'same line=solid,2 cpts magenta nsta notit'

endscript

defscript pldishes 

  Dx = 50
  Dy = 50
  calcul 
  plot2d pknD50m  k/$Da pnoise*$PNOISE/100 fracmodok>0.2 'same line=solid,2 cpts skyblue nsta notit'

  Dx = 75
  Dy = 75
  calcul 
  plot2d pknD75m  k/$Da pnoise*$PNOISE/100 fracmodok>0.2 'same line=solid,2 cpts blue nsta notit'

  Dx = 100
  Dy = 100
  calcul 
  plot2d pknD100m  k/$Da pnoise*$PNOISE/100 fracmodok>0.2 'same line=solid,2 cpts navyblue nsta notit'

  Dx = 200
  Dy = 200
  calcul 
  plot2d pknD200m  k/$Da pnoise*$PNOISE/100 fracmodok>0.2 'same line=solid,2 cpts blueviolet nsta notit'

  set lines (  'D:Dish50m' 'D:Dish75m' 'D:Dish100m' 'D:Dish200m' )
  set cols ( skyblue blue navyblue blueviolet )
  textdrawer lines cols 'frame font=helvetica,bold,16  inset=0.1,0.3,0.15,0.35' 
  settitle ' PNoise(k) : Dish D=50m,75m,100m,200m' ' ' 'font=helvetica,bold,16' 
endscript



defscript plnoisenancay
#  Dx = 5*0.95
#  Dy = 5*0.95
#  calcul 
#  plot2d pknnan24 x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts red nsta notit'
#  plot2d pknnan25 x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts orange nsta notit'
#  Dx = 3.5*0.95
#  Dy = 3.5*0.95
#  calcul 
#  plot2d pknnan36 x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts red nsta notit'
#  plot2d pknnan40 x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts violet nsta notit'
#  Dx = 7*0.9
#  Dy = 2*$Lambda*0.9
#  calcul 
#  plot2d pknpit2cyl x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts green nsta notit'
  Dx = 9*0.9
  Dy = 2*$Lambda*0.9
  calcul 
  plot2d pknpit2cylw x/$Da val*$PNOISE nb>10 'same line=solid,2 cpts red nsta notit'
  set lines ( '2 Cylinders 20mx9m' )
  set cols ( red )
  textdrawer lines cols 'noframe font=helvetica,bold,16 inset=0.4,0.4,0.4,0.4' 

endscript

defscript ABCD
#  plnoisedish
#  zone 1 2
#  newwin 1 1 1000 600
  plpklss
  plnoiseABCD
#  plfracABCD
endscript 

defscript EFGH
#  plnoisedish
#  zone 1 2
#  newwin 1 1 900 600
  plpklss
  plnoiseEFGH
#  plfracEFGH
endscript 

defscript DDDD
#  plnoisedish
#  zone 1 2
  plpklss
  plnoiseDishes
#  plfracABCD
endscript 

defscript BB
  plnoisedish
  plpklss
  plnoiseB
endscript 

defscript CC
  plnoisedish
  plpklss
  plnoiseC2
endscript 

defscript nancay
  plnoisedish
  plpklss
  plnoisenancay
endscript 

defscript POSCOV
  openppf hdt_repnan128.ppf
  rename mdish mdB128d
  rename h2rep uvB128d

  openppf hdt_repf11x11.ppf
  rename mdish mdA121d
  rename h2rep uvA121d

  openppf hdt_repconfC.ppf
  rename mdish mdC129d
  rename h2rep uvC129d

  newwin 2 1 800 400
  zone 2  1
  nt2d mdB128d posx posy - - - - 'xylimits=-10,90,-10,90 marker=circle,15 notit nsta red '
  setaxelabels ' X (meters) ' ' Y (meters) ' 'font=helvetica,bolditalic,16' 
  settitle '(b) 128 D=5m dishes in 8 rows ' ' ' 'font=helvetica,bold,16' 
  nt2d mdC129d posx posy - - - - 'xylimits=-10,90,-10,90 marker=circle,15 notit nsta red'
  setaxelabels ' X (meters) ' ' Y (meters) ' 'font=helvetica,bolditalic,16' 
  settitle '(c) 129 D=5m dishes ' ' ' 'font=helvetica,bold,16' 
  pssetfilename configab.ps
  w2ps
  w2eps configab.eps
  psclosefile 

  newwin 2 2 800 800 
  disp uvA121d 'h2disp=img colbr128 h2dyn=1,80 notit nsta'
  setaxelabels 'u (Radian^-1) u=1000->~21 arcmin' ' v (Radian^-1) ' 'font=helvetica,bolditalic,12' 
  settitle 'u-v coverage, (a) 11x11 D=5m dishes Over 55mx55m' ' ' 'font=helvetica,bold,12' 
  disp uvB128d 'h2disp=img colbr128 h2dyn=1,80 notit nsta'
  setaxelabels 'u (Radian^-1) u=1000->~21 arcmin' ' v (Radian^-1) ' 'font=helvetica,bolditalic,12' 
  settitle 'u-v coverage, (b) 8 row of 16xD=5m dishes Over 80mx80m' ' ' 'font=helvetica,bold,12' 
  disp uvC129d 'h2disp=img colbr128 h2dyn=1,80 notit nsta'
  setaxelabels 'u (Radian^-1) u=1000->~21 arcmin' ' v (Radian^-1) ' 'font=helvetica,bolditalic,12' 
  settitle 'u-v coverage, (c) 129 D=5m dishes Over 80mx80m' ' ' 'font=helvetica,bold,12' 
  disp h2drD75m 'h2disp=img colbr128 h2dyn=0.02,0.8 notit nsta'
  setaxelabels 'u (Radian^-1) u=1000->~21 arcmin' ' v (Radian^-1) ' 'font=helvetica,bolditalic,12' 
  settitle 'u-v coverage, (d) D=75 m dish with 100 beams' ' ' 'font=helvetica,bold,12' 

  pssetfilename uvcovabcd.ps
  w2ps
  w2eps uvcovabcd.eps
  psclosefile 

#  imag uvB128d 'colbr128 stdaxes showcmap=right' 
  
endscript