-
Notifications
You must be signed in to change notification settings - Fork 5
/
kinetic_h2.pro
2263 lines (2230 loc) · 93.5 KB
/
kinetic_h2.pro
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
;+
; Kinetic_H2.pro
;
; This subroutine is part of the "KN1D" atomic and molecular neutal transport code.
;
; This subroutine solves a 1-D spatial, 2-D velocity kinetic neutral transport
; problem for molecular hydrogen or deuterium (H2) by computing successive generations of
; charge exchange and elastic scattered neutrals. The routine handles electron-impact
; ionization and dissociation, molecular ion charge exchange, and elastic
; collisions with hydrogenic ions, neutral atoms, and molecules.
;
; The positive vx half of the atomic neutral distribution function is inputted at x(0)
; (with arbitrary normalization). The desired flux on molecules entering the slab geometry at
; x(0) is specified. Background profiles of plasma ions, (e.g., Ti(x), Te(x), n(x), vxi(x),...)
; and atomic distribution function (fH) is inputted. (fH can be computed by procedure
; "Kinetic_H.pro".) Optionally, a molecular source profile (SH2(x)) is also inputted.
; The code returns the molecular hydrogen distribution function, fH2(vr,vx,x) for all
; vx, vr, and x of the specified vr,vx,x grid. The atomic (H) and ionic (P) hydrogen
; source profiles and the atomic source velocity distribution functions
; resulting from Franck-Condon reaction product energies of H are also returned.
;
; Since the problem involves only the x spatial dimension, all distribution functions
; are assumed to have rotational symmetry about the vx axis. Consequently, the distributions
; only depend on x, vx and vr where vr =sqrt(vy^2+vz^2)
;
; History:
;
; B. LaBombard First coding based on Kinetic_Neutrals.pro 22-Dec-2000
;
; For more information, see write-up: "A 1-D Space, 2-D Velocity, Kinetic
; Neutral Transport Algorithm for Hydrogen Molecules in an Ionizing Plasma", B. LaBombard
;
; Note: Variable names contain characters to help designate species -
; atomic neutral (H), molecular neutral (H2), molecular ion (HP), proton (i) or (P)
;
;________________________________________________________________________________
pro Kinetic_H2,vx,vr,x,Tnorm,mu,Ti,Te,n,vxi,fH2BC,GammaxH2BC,NuLoss,PipeDia,fH,SH2,$
fH2,nH2,GammaxH2,VxH2,pH2,TH2,qxH2,qxH2_total,Sloss,QH2,RxH2,QH2_total,AlbedoH2,WallH2,$
truncate=truncate,$
nHP,THP,fSH,SH,SP,SHP,NuE,NuDis,$
Simple_CX=Simple_CX,Max_Gen=Max_Gen,Compute_H_Source=Compute_H_Source,$
No_Sawada=No_Sawada,H2_H2_EL=H2_H2_EL,H2_P_EL=H2_P_EL,H2_H_EL=_H2_H_EL,H2_HP_CX=H2_HP_CX,ni_correct=ni_correct,$
ESH=ESH,Eaxis=Eaxis,error=error,compute_errors=compute_errors,$
plot=plot,debug=debug,debrief=debrief,pause=pause
common Kinetic_H2_Output,piH2_xx,piH2_yy,piH2_zz,RxH2CX,RxH_H2,RxP_H2,RxW_H2,EH2CX,EH_H2,EP_H2,EW_H2,Epara_PerpH2_H2
common Kinetic_H2_Errors,Max_dx,vbar_error,mesh_error,moment_error,C_Error,CX_Error,Swall_Error,H2_H2_error,Source_Error,$
qxH2_total_error,QH2_total_error
;
; Input:
; vx(*) - fltarr(nvx), normalized x velocity coordinate
; [negative values, positive values],
; monotonically increasing. Note: a nonuniform mesh can be used.
; Dimensional velocity (note: Vth is based on ATOM mass)
; is v = Vth * vx where Vth=sqrt(2 k Tnorm/(mH*mu))
; Note: nvx must be even and vx(*) symmetric about
; zero but not contain a zero element
; vr(*) - fltarr(nvr), normalized radial velocity coordinate
; [positive values], monotonically increasing. Note: a non-uniform mesh can be used.
; Dimensional velocity is v = Vth * vr where Vth=sqrt(2 k Tnorm/(mH*mu))
; Note: vr must not contain a zero element
; x(*) - fltarr(nx), spatial coordinate (meters),
; positive, monontonically increasing. Note: a non-uniform mesh can be used.
; Tnorm - Float, temperature corresponding to the thermal speed (see vx and vr above) (eV)
; mu - Float, 1=hydrogen, 2=deuterium
; Ti - fltarr(nx), Ion temperature profile (eV)
; Te - fltarr(nx), electron temperature profile (eV)
; n - fltarr(nx), electron density profile (m^-3)
; vxi - fltarr(nx), x-directed plasma ion and molecular ion flow profile (m s^-1)
; fH2BC - fltarr(nvr,nvx), this is an input boundary condition
; specifying the shape of the neutral molecule velocity distribution
; function at location x(0). Normalization is arbitrary.
; Only values with positive vx, fH2BC(*,nvx/2:*) are used
; by the code.
; GammaxH2BC - float, desired neutral molecule flux density in the +Vx
; direction at location x(0) (m^-2 s^-1)
; fH2BC is scaled to yield this flux density.
; NuLoss - fltarr(nx), characteristic loss frequency for HP molecules (1/s)
; (for open field lines, this is ~Cs/L). If this variable is undefined,
; then NuLoss set set to zero.
; PipeDia - fltarr(nx), effective pipe diameter (meters)
; This variable allows collisions with the 'side-walls' to be simulated.
; If this variable is undefined, then PipeDia set set to zero. Zero values
; of PipeDia are ignored (i.e., treated as an infinite diameter).
; fH - fltarr(nvr,nvx,nx), neutral atomic velocity distribution
; function. fH is normalized so that the atomic neutral density, nH(k), is
; defined as the velocity space integration:
; nH(k)=total(Vr2pidVr*(fH(*,*,k)#dVx))
; If this variable is undefined, then no molecule-atom collisions are included.
; NOTE: dVx is velocity space differential for Vx axis and Vr2pidVr = Vr*!pi*dVr
; with dVr being velocity space differential for Vr axis.
; SH2 - fltarr(nx), source profile of wall-temperature (room temperature) H2 molecules (m^-3 s^-1)
; If this variable is undefined, it is set equal to zero.
;
; Input & Output:
; fH2 - fltarr(nvr,nvx,nx), neutral molecule velocity distribution
; function.
; 'Seed' values for this may be specified on input. If this parameter
; is undefined, then a zero 'seed' value will be used.
; The algorithm outputs a self-consistent fH2.
; fH2 is normalized so that the neutral density, nH2(k), is defined as
; the velocity space integration:
; nH2(k)=total(Vr2pidVr*(fH2(*,*,k)#dVx))
; nHP - fltarr(nx), molecular ion density profile (m^-3)
; 'Seed' values for this may be specified on input. If this parameter
; is undefined, then a zero 'seed' value will be used.
; The algorithm outputs a self-consistent profile for nHP.
; THP - fltarr(nx), molecular ion temperature profile (m^-3)
; 'Seed' values for this may be specified on input. If this parameter
; is undefined, then a 'seed' value of 3 eV will be used.
; The algorithm outputs a self-consistent profile for THP.
;
; Output:
; fH2 - fltarr(nvr,nvx,nx), neutral molecule velocity distribution
; function. fH2 is normalized so that the neutral density, nH2(k), is
; defined as the velocity space integration:
; nH2(k)=total(Vr2pidVr*(fH2(*,*,k)#dVx))
; nH2 - fltarr(nx), neutral molecule density profile (m^-3)
; GammaxH2 - fltarr(nx), neutral flux profile (# m^-2 s^-1)
; computed from GammaxH2(k)=Vth*total(Vr2pidVr*(fH2(*,*,k)#(Vx*dVx)))
; VxH2 - fltarr(nx), neutral velocity profile (m s^-1)
; computed from GammaxH2/nH2
;
; To aid in computing the some of the quantities below, the procedure internally
; defines the quantities:
; vr2vx2_ran(i,j,k)=vr(i)^2+(vx(j)-VxH2(k))^2
; which is the magnitude of 'random v^2' at each mesh point
; vr2vx2(i,j,k)=vr(i)^2+vx(j)^2
; which is the magnitude of 'total v^2' at each mesh point
; q=1.602177D-19, mH=1.6726231D-27
; C(*,*,*) is the right hand side of the Boltzmann equation, evaluated
; using the computed neutral distribution function
;
; pH2 - fltarr(nx), neutral pressure (eV m^-2) computed from:
; pH2(k)~vth2*total(Vr2pidVr*(vr2vx2_ran(*,*,k)*fH2(*,*,k))#dVx))*(mu*mH)/(3*q)
; TH2 - fltarr(nx), neutral temperature profile (eV) computed from: TH2=pH2/nH2
; qxH2 - fltarr(nx), neutral random heat flux profile (watts m^-2) computed from:
; qxH2(k)~vth3*total(Vr2pidVr*((vr2vx2_ran(*,*,k)*fH2(*,*,k))#(dVx*(vx-VxH2(k)))))*0.5*(mu*mH)
; qxH2_total - fltarr(nx), total neutral heat flux profile (watts m^-2)
; This is the total heat flux transported by the neutrals:
; qxH2_total=(0.5*nH2*(mu*mH)*VxH2*VxH2 + 2.5*pH2*q)*VxH2 + piH2_xx*VxH2 + qxH2
; Sloss - fltarr(nx), H2 loss rate from ionization and dissociation (SH2) (m^-3 s^-1)
; QH2 - fltarr(nx), rate of net thermal energy transfer into neutral molecules (watts m^-3) computed from
; QH2(k)~vth2*total(Vr2pidVr*((vr2vx2_ran(*,*,k)*C(*,*,k))#dVx))*0.5*(mu*mH)
; RxH2 - fltarr(nx), rate of x momentum transfer to neutral molecules (=force, N m^-2).
; RxH2(k)~Vth*total(Vr2pidVr*(C(*,*,k)#(dVx*(vx-VxH2(k)))))*(mu*mH)
; QH2_total - fltarr(nx), net rate of total energy transfer into neutral molecules
; = QH2 + RxH2*VxH2 - 0.5*(mu*mH)*(Sloss-SH2)*VxH2*VxH2 (watts m^-3)
; AlbedoH2 - float, Ratio of molceular particle flux with Vx < 0 divided by particle flux
; with Vx > 0 at x=x(0)
; (Note: For SH2 non-zero, the flux with Vx < 0 will include
; contributions from molecular hydrogen sources within the 'slab'.
; In this case, this parameter does not return the true 'Albedo'.)
; WallH2 - fltarr(nx), molecular sink rate and source rate from interation with 'side walls' (m^-3 s^-1)
; Note: There is no loss of molecules when they strike the side walls in this model. The
; molecules are just relaunched with a Maxwellian distribution at the wall temperature.
;
; nHP - fltarr(nx), molecular ion density profile (m^-3)
; 'Seed' values for this may be specified on input.
; THP - fltarr(nx), molecular ion temperature profile (m^-3)
; 'Seed' values for this may be specified on input.
; fSH - fltarr(nvr,nvx,nx), H source velocity distribution.
; fSH is normalized so that the total atomic neutral
; source, SH(k), is defined as the velocity space integration:
; SH(k)=total(Vr2pidVr*(fSH(*,*,k)#dVx))
; SH - fltarr(nx), atomic neutral source profile (m^-3 s^-1)
; computed from total(Vr2pidVr*(fSH(*,*,k)#dVx))
; SP - fltarr(nx), proton source profile (m^-3 s^-1)
; SHP - fltarr(nx), molecular ion source profile (m^-3 s^-1)
; NuE - fltarr(nx), energy equilibration frequency of molecular ion with bulk plasma (1/s)
; NuDis - fltarr(nx), molecular ion dissociation frequency (1/s)
;
; KEYWORDS:
; Output:
; ESH - fltarr(nvr,nx) - returns normalized H source energy distribution
; derived from velocity-space distribution:
; ESH(*,*) = 0.5*mu*mH*Vr^2 FSH(*,vx=0,*)*4*pi*Vr/(mu*mH)
; Eaxis - fltarr(nvr) - returns energy coordinate for ESH (eV) = 0.5*mu*mH*Vr^2
; error - Returns error status: 0=no error, solution returned
; 1=error, no solution returned
;
; COMMON BLOCK KINETIC_H2_OUTPUT
; Output:
; piH2_xx - fltarr(nx), xx element of stress tensor (eV m^-2) computed from:
; piH2_xx(k)~vth2*total(Vr2pidVr*(fH2(*,*,k)#(dVx*(vx-VxH2(k))^2)))*(mu*mH)/q - pH2
; piH2_yy - fltarr(nx), yy element of stress tensor (eV m^-2) computed from:
; piH2_yy(k)~vth2*total((Vr2pidVr*Vr^2)*(fH2(*,*,k)#dVx))*(mu*mH)/q - pH2
; piH2_zz - fltarr(nx), zz element of stress tensor (eV m^-2) = piH2_yy
; Note: cylindrical system relates r^2 = y^2 + z^2. All other stress tensor elements are zero.
;
; The following momentum and energy transfer rates are computed from charge-exchange collsions between species:
; RxH2CX - fltarr(nx), rate of x momentum transfer from molecular ions to neutral molecules (=force/vol, N m^-3).
; EH2CX - fltarr(nx), rate of energy transfer from molecular ions to neutral molecules (watts m^-3).
;
; The following momentum and energy transfer rates are computed from elastic collsions between species:
; RxH_H2 - fltarr(nx), rate of x momentum transfer from neutral atoms to molecules (=force/vol, N m^-3).
; RxP_H2 - fltarr(nx), rate of x momentum transfer from hydrogen ions to neutral molecules (=force/vol, N m^-3).
; EH_H2 - fltarr(nx), rate of energy transfer from neutral atoms to molecules (watts m^-3).
; EP_H2 - fltarr(nx), rate of energy transfer from hydrogen ions to neutral molecules (watts m^-3).
;
; The following momentum and energy transfer rates are computed from collisions with the 'side-walls'
; RxW_H2 - fltarr(nx), rate of x momentum transfer from wall to molecules (=force/vol, N m^-3).
; EW_H2 - fltarr(nx), rate of energy transfer from wall to neutral molecules (watts m^-3).
;
; The following is the rate of parallel to perpendicular energy transfer computed from elastic collisions
; Epara_PerpH2_H2 - fltarr(nx), rate of parallel to perp energy transfer within molecular hydrogen species (watts m^-3).
;
; KEYWORDS:
; Input:
; truncate - float, stop computation when the maximum
; increment of neutral density normalized to
; inputed neutral density is less than this
; value in a subsequent generation. Default value is 1.0e-4
;
; Simple_CX - if set, then use CX source option (B): Neutral molecules are born
; in velocity with a distribution proportional to the local
; molecular ion distribution function. Simple_CX=1 is default.
;
; if not set, then use CX source option (A): The CX source
; neutral molecule distribution function is computed by evaluating the
; the CX cross section for each combination of (vr,vx,vr',vx')
; and convolving it with the molecule neutral distribution function.
; This option requires more CPU time and memory.
;
; Max_gen - integer, maximum number of collision generations to try including before giving up.
; Default is 50.
; Compute_H_Source - if set, then compute fSH, SH, SP, and SHP
;
; No_Sawada - if set, then DO NOT correct reaction rates according to
; results from collisional-radiative model of Sawada
; [Sawada, K. and Fujimoto, T., Journal of Applied Physics 78 (1995) 2913.]
; H2_H2_EL - if set, then include H2 -> H2 elastic self collisions
; Note: if H2_H2_EL is set, then algorithm iterates fH2 until
; self consistent fH2 is achieved.
; H2_HP_CX - if set, then include H2 -> H2(+) charge exchange collisions
; Note: if H2_HP_CX is set, then algorithm iterates until
; self consistent nHp is achieved.
; H2_P_EL - if set, then include H2 -> H(+) elastic collisions (does not lead to an iteration loop)
; H2_H_EL - if set, then include H2 -> H elastic collisions (does not lead to an iteration loop)
; ni_correct - if set, then algorithm corrects hydrogen ion density
; according to quasineutrality: ni=ne-nHp
; This is done in an iteration loop.
;
; Compute_Errors - if set, then return error estimates in common block KINETIC_H2_ERRORS below
;
; plot - 0= no plots, 1=summary plots, 2=detail plots, 3=very detailed plots
; debug - 0= do not execute debug code, 1=summary debug, 2=detail debug, 3=very detailed debug
; debrief - 0= do not print, 1=print summary information, 2=print detailed information
; pause - if set, then pause between plots
;
; COMMON BLOCK KINETIC_H2_ERRORS
;
; if COMPUTE_ERRORS keyword is set then the following is returned in common block KINETIC_H2_ERRORS
;
; Max_dx - float(nx), Max_dx(k) for k=0:nx-2 returns maximum
; allowed x(k+1)-x(k) that avoids unphysical negative
; contributions to fH2
; Vbar_error - float(nx), returns numerical error in computing
; the speed of neutrals averged over maxwellian distribution;
; over a temperature range spanning the Franck Condon energies
; of reactions R2, R3, R4, R5, R6, R7, R8, R10
; The average speed should be:
; vbar_exact=2*Vth*sqrt(Ti(*)/Tnorm)/sqrt(!pi)
; Vbar_error returns: abs(vbar-vbar_exact)/vbar_exact
; where vbar is the numerically computed value.
; mesh_error - fltarr(nvr,nvx,nx), normalized error of solution
; based on substitution into Boltzmann equation.
; moment_error - fltarr(nx,m), normalized error of solution
; based on substitution into velocity space
; moments (v^m) of Boltzmann equation, m=[0,1,2,3,4]
; C_error - fltarr(nx), normalized error in charge exchange and elastic scattering collision
; operator. This is a measure of how well the charge exchange and
; elastic scattering portions of the collision operator
; conserve particles.
; CX_error - fltarr(nx), normalized particle conservation error in charge exchange collision operator.
; H2_H2_error - fltarr(nx,[0,1,2]) return normalized errors associated with
; particle [0], x-momentum [1], and total energy [2] convervation of the elastic self-collision operator
;
; Source_error - fltarr(nx), normalized error in mass continuity equation. This is a measure
; of how well atomic plus molecular mass continuity is satisfied.
; qxH2_total_error - fltarr(nx), normalized error estimate in computation of qxH2_total
; QH2_total_error - fltarr(nx), normalized error estimate in computation of QH2_total
;
; History:
; 22-Dec-2000 - B. LaBombard - first coding.
; 06-Feb-2001 - B. LaBombard - added elastic collisions and molecular sources
;
;______________________________________________________________________
;-
prompt='Kinetic_H2 => '
;
common Kinetic_H2_input,vx_s,vr_s,x_s,Tnorm_s,mu_s,Ti_s,Te_s,n_s,vxi_s,fH2BC_s,GammaxH2BC_s,NuLoss_s,PipeDia_s,fH_s,SH2_s,fH2_s,$
nHP_s,THP_s,Simple_CX_s,Sawada_s,H2_H2_EL_s,H2_P_EL_s,H2_H_EL_s,H2_HP_CX_s,ni_correct_s
common Kinetic_H2_internal,vr2vx2,vr2vx_vxi2,fw_hat,fi_hat,fHp_hat,EH2_P,sigv,Alpha_Loss,v_v2,v_v,vr2_vx2,vx_vx,$
Vr2pidVrdVx,SIG_CX,SIG_H2_H2,SIG_H2_H,SIG_H2_P,Alpha_CX,Alpha_H2_H,MH2_H2_sum,Delta_nH2s
common Kinetic_H2_H_moments,nH,VxH,TH
;
; Internal Debug switches
;
shifted_Maxwellian_debug=0
CI_Test=1
Do_Alpha_CX_Test=0
;
; Internal Tolerances
;
DeltaVx_tol=.01
Wpp_tol=.001
;
; Test input parameters
;
key_default,truncate,1.0e-4
key_default,Compute_Errors,0
key_default,Compute_H_Source,0
key_default,plot,0
key_default,debug,0
key_default,pause,0
key_default,debrief,0
if debug gt 0 then plot=plot > 1
if debug gt 0 then debrief=debrief > 1
if debug gt 0 then pause=1
key_default,Simple_CX,1
key_default,Max_Gen,50
key_default,No_Sawada,0
if No_Sawada eq 0 then Sawada=1 else Sawada=0
key_default,H2_H2_EL,0
key_default,H2_P_EL,0
key_default,_H2_H_EL,0
key_default,H2_HP_CX,0
key_default,ni_correct,0
error=0
nvr=n_elements(vr)
nvx=n_elements(vx)
nx=n_elements(x)
vr=double(vr)
vx=double(vx)
x=double(x)
dx=x-shift(x,1) & dx=dx(1:*)
notpos=where(dx le 0.0,count)
if count gt 0 then begin & print,prompt+'x(*) must be increasing with index!' & error=1 & goto,return & endif
if (nvx mod 2) ne 0 then begin & print,prompt+'Number of elements in vx must be even!' & error=1 & goto,return & endif
if n_elements(Ti) ne nx then begin & print,prompt+'Number of elements in Ti and x do not agree!' & error=1 & goto,return & endif
if type_of(vxi) eq 0 then vxi=dblarr(nx)
if n_elements(vxi) ne nx then begin & print,prompt+'Number of elements in vxi and x do not agree!' & error=1 & goto,return & endif
if n_elements(Te) ne nx then begin & print,prompt+'Number of elements in Te and x do not agree!' & error=1 & goto,return & endif
if n_elements(n) ne nx then begin & print,prompt+'Number of elements in n and x do not agree!' & error=1 & goto,return & endif
if type_of(NuLoss) eq 0 then NuLoss=fltarr(nx)
if n_elements(NuLoss) ne nx then begin & print,prompt+'Number of elements in NuLoss and x do not agree!' & error=1 & goto,return & endif
if type_of(PipeDia) eq 0 then PipeDia=dblarr(nx)
if n_elements(PipeDia) ne nx then begin & print,prompt+'Number of elements in PipeDia and x do not agree!' & error=1 & goto,return & endif
if type_of(GammaxH2BC) eq 0 then begin & print,prompt+'GammaxH2BC is not defined!' & error=1 & goto,return & endif
if type_of(fH) eq 0 then fH=dblarr(nvr,nvx,nx)
if n_elements(fH(*,0,0)) ne nvr then begin & print,prompt+'Number of elements in fH(*,0,0) and vr do not agree!' & error=1 & goto,return & endif
if n_elements(fH(0,*,0)) ne nvx then begin & print,prompt+'Number of elements in fH(0,*,0) and vx do not agree!' & error=1 & goto,return & endif
if n_elements(fH(0,0,*)) ne nx then begin & print,prompt+'Number of elements in fH(0,0,*) and x do not agree!' & error=1 & goto,return & endif
if n_elements(fH2BC(*,0)) ne nvr then begin & print,prompt+'Number of elements in fH2BC(*,0) and vr do not agree!' & error=1 & goto,return & endif
if n_elements(fH2BC(0,*)) ne nvx then begin & print,prompt+'Number of elements in fH2BC(0,*) and vx do not agree!' & error=1 & goto,return & endif
if type_of(fH2) eq 0 then fH2=dblarr(nvr,nvx,nx)
if n_elements(fH2(*,0,0)) ne nvr then begin & print,prompt+'Number of elements in fH2(*,0,0) and vr do not agree!' & error=1 & goto,return & endif
if n_elements(fH2(0,*,0)) ne nvx then begin & print,prompt+'Number of elements in fH2(0,*,0) and vx do not agree!' & error=1 & goto,return & endif
if n_elements(fH2(0,0,*)) ne nx then begin & print,prompt+'Number of elements in fH2(0,0,*) and x do not agree!' & error=1 & goto,return & endif
if type_of(SH2) eq 0 then SH2=dblarr(nx)
if n_elements(SH2) ne nx then begin & print,prompt+'Number of elements in SH2 and x do not agree!' & error=1 & goto,return & endif
if type_of(nHP) eq 0 then nHP=dblarr(nx)
if n_elements(nHP) ne nx then begin & print,prompt+'Number of elements in nHP and x do not agree!' & error=1 & goto,return & endif
if type_of(THP) eq 0 then THP=dblarr(nx)+3.0
if n_elements(THP) ne nx then begin & print,prompt+'Number of elements in THP and x do not agree!' & error=1 & goto,return & endif
if total(abs(vr)) le 0.0 then begin & print,prompt+'vr is all 0!' & error=1 & goto,return & endif
ii=where(vr le 0,count)
if count gt 0 then begin & print,prompt+'vr contains zero or negative element(s)!' & error=1 & goto,return & endif
if total(abs(vx)) le 0.0 then begin & print,prompt+'vx is all 0!' & error=1 & goto,return & endif
if total(x) le 0.0 then begin & print,prompt+'Total(x) is less than or equal to 0!' & error=1 & goto,return & endif
if type_of(Tnorm) eq 0 then begin & print,prompt+'Tnorm is not defined!' & error=1 & goto,return & endif
if type_of(mu) eq 0 then begin & print,prompt+'mu is not defined!' & error=1 & goto,return & endif
if mu ne 1 and mu ne 2 then begin & print,prompt+'mu must be 1 or 2!' & error=1 & goto,return & endif
_e='e!U-!N'
if mu eq 1 then begin
_p='H!U+!N'
_H='H!U0!N'
_H1s='H(1s)'
_H2s='H!U*!N(2s)'
_H2p='H!U*!N(2p)'
_Hn2='H!U*!N(n=2)'
_Hn3='H!U*!N(n=3)'
_Hn='H!U*!N(n>=2)'
_HH='H!D2!N'
_Hp='H!D2!U+!N'
endif else begin
_p='D!U+!N'
_H='D!U0!N'
_H1s='D(1s)'
_H2s='D!U*!N(2s)'
_H2p='D!U*!N(2p)'
_Hn2='D!U*!N(n=2)'
_Hn3='D!U*!N(n=3)'
_Hn='D!U*!N(n>=2)'
_HH='D!D2!N'
_Hp='D!D2!U+!N'
endelse
plus=' + '
arrow=' -> '
elastic=' (elastic)'
_R1=_e+plus+_HH+arrow+_e+plus+_Hp+plus+_e
_R2=_e+plus+_HH+arrow+_e+plus+_H1s+plus+_H1s
_R3=_e+plus+_HH+arrow+_e+plus+_H1s+plus+_H2s
_R4=_e+plus+_HH+arrow+_e+plus+_p+plus+_H1s+plus+_e
_R5=_e+plus+_HH+arrow+_e+plus+_H2p+plus+_H2s
_R6=_e+plus+_HH+arrow+_e+plus+_H1s+plus+_Hn3
_R7=_e+plus+_Hp+arrow+_e+plus+_p+plus+_H1s
_R8=_e+plus+_Hp+arrow+_e+plus+_p+plus+_Hn2
_R9=_e+plus+_Hp+arrow+_e+plus+_p+plus+_p+plus+_e
_R10=_e+plus+_Hp+arrow+_H1s+plus+_Hn
_R11=_HH+plus+_p+arrow+_HH+plus+_p+elastic
_R12=_HH+plus+_H+arrow+_HH+plus+_H+elastic
_R13=_HH+plus+_HH+arrow+_HH+plus+_HH+elastic
_R14=_HH+plus+_Hp+arrow+_Hp+plus+_HH
_Rn=[' ',_R1,_R2,_R3,_R4,_R5,_R6,_R7,_R8,_R9,_R10,_R11,_R12,_R13,_R14]
in=where(vx lt 0,count)
if count lt 1 then begin & print,prompt+'vx contains no negative elements!' & error=1 & goto,return & endif
ip=where(vx gt 0,count)
if count lt 1 then begin & print,prompt+'vx contains no positive elements!' & error=1 & goto,return & endif
iz=where(vx eq 0,count)
if count gt 0 then begin & print,prompt+'vx contains one or more zero elements!' & error=1 & goto,return & endif
diff=where(vx(ip) ne -reverse(vx(in)),count)
if count gt 0 then begin & print,prompt+'vx array elements are not symmetric about zero!' & error=1 & goto,return & endif
fH2BC_input=fH2BC
fH2BC_input(*)=0.0
fH2BC_input(*,ip)=fH2BC(*,ip)
test=total(fH2BC_input)
if test le 0.0 then begin & print,prompt+'Values for fH2BC(*,*) with vx > 0 are all zero!' & error=1 & goto,return & endif
;
; Output variables
;
nH2=dblarr(nx)
GammaxH2=dblarr(nx)
VxH2=dblarr(nx)
pH2=dblarr(nx)
TH2=dblarr(nx)
NuDis=dblarr(nx)
NuE=dblarr(nx)
qxH2=dblarr(nx)
qxH2_total=dblarr(nx)
Sloss=dblarr(nx)
WallH2=dblarr(nx)
QH2=dblarr(nx)
RxH2=dblarr(nx)
QH2_total=dblarr(nx)
piH2_xx=dblarr(nx)
piH2_yy=dblarr(nx)
piH2_zz=dblarr(nx)
RxH2CX=dblarr(nx)
RxH_H2=dblarr(nx)
RxP_H2=dblarr(nx)
RxW_H2=dblarr(nx)
EH2CX=dblarr(nx)
EH_H2=dblarr(nx)
EP_H2=dblarr(nx)
EW_H2=dblarr(nx)
Epara_PerpH2_H2=dblarr(nx)
AlbedoH2=0.0D0
fSH=dblarr(nvr,nvx,nx)
SH=dblarr(nx)
SP=dblarr(nx)
SHP=dblarr(nx)
ESH=dblarr(nvr,nx)
Eaxis=dblarr(nx)
;
; Internal variables
;
mH=1.6726231D-27
q=1.602177D-19
k_boltz=1.380658D-23 &;Bolzmann's constant, J K^-1
Twall=293.0*k_boltz/q &;room temperature (eV)
Work=dblarr(nvr*nvx)
fH2G=dblarr(nvr,nvx,nx)
NH2G=dblarr(nx,max_gen+1)
Vth=sqrt(2*q*Tnorm/(mu*mH))
Vth2=vth*vth
Vth3=Vth2*Vth
fH2s=dblarr(nx)
nH2s=dblarr(nx)
THPs=dblarr(nx)
nHPs=dblarr(nx)
Alpha_H2_H2=dblarr(nvr,nvx)
Omega_H2_P=dblarr(nx)
Omega_H2_H=dblarr(nx)
Omega_H2_H2=dblarr(nx)
VxH2G=dblarr(nx)
TH2G=dblarr(nx)
Wperp_paraH2=dblarr(nx)
vr2vx2_ran2=dblarr(nvr,nvx)
vr2_2vx_ran2=dblarr(nvr,nvx)
vr2_2vx2_2D=dblarr(nvr,nvx)
RxCI_CX=dblarr(nx)
RxCI_H_H2=dblarr(nx)
RxCI_P_H2=dblarr(nx)
Epara_Perp_CI=dblarr(nx)
CI_CX_error=fltarr(nx)
CI_H_H2_error=fltarr(nx)
CI_P_H2_error=fltarr(nx)
CI_H2_H2_error=fltarr(nx)
Maxwell=dblarr(nvr,nvx,nx)
Make_dVr_dVx,vr,vx,Vr2pidVr,VrVr4pidVr,dVx,vrL=vrL,vrR=vrR,vxL=vxL,vxR=vxR,$
Vol=Vol,Vth_DeltaVx=Vth_DVx,Vx_DeltaVx=Vx_DVx,Vr_DeltaVr=Vr_DVr,Vr2Vx2=Vr2Vx2_2D,$
jpa=jpa,jpb=jpb,jna=jna,jnb=jnb
;
; Vr^2-2*Vx^2
;
for i=0,nvr-1 do vr2_2vx2_2D(i,*)=vr(i)^2-2*vx^2
;
; Theta-prime coordinate
;
ntheta=5 &; use 5 theta mesh points for theta integration
dTheta=replicate(1.0d0,ntheta)/ntheta
theta=!pi*(dindgen(ntheta)/ntheta+0.5/ntheta)
cos_theta=cos(theta)
;
; Determine energy space differentials
;
Eaxis=vth2*0.5*mu*mH*vr^2/q
_Eaxis=[Eaxis,2*Eaxis(nvr-1)-Eaxis(nvr-2)]
Eaxis_mid=[0.0,0.5*(_Eaxis+shift(_Eaxis,-1))]
dEaxis=shift(Eaxis_mid,-1)-Eaxis_mid
dEaxis=dEaxis(0:nvr-1)
;
; Scale input molecular distribution function to agree with desired flux
;
gamma_input=1.0
if abs(GammaxH2BC) gt 0 then gamma_input=vth*total(Vr2pidVr*(fH2BC_input#(Vx*dVx)))
ratio=abs(GammaxH2BC)/gamma_input
fH2BC_input=fH2BC_input*ratio
if abs(ratio-1) gt 0.01*truncate then begin
fH2BC=fH2BC_input
endif
fH2(*,ip,0)=fH2BC_input(*,ip)
;
; if fH is zero, then turn off elastic H2 <-> H collisions
;
H2_H_EL=_H2_H_EL
if total(fH) le 0.0 then H2_H_EL=0
;
; Set iteration scheme
;
fH2_iterate=0
if (H2_H2_EL ne 0) or (H2_HP_CX ne 0) or (H2_H_EL ne 0) or (H2_P_EL ne 0) or (ni_correct ne 0) then fH2_iterate=1
fH2_generations=0
if (fH2_iterate ne 0) then fH2_generations=1
;
; Set flags to make use of previously computed local parameters
;
New_Grid=1
if type_of(vx_s) ne 0 then begin
test=0
ii=where(vx_s ne vx,count) & test=test+count
ii=where(vr_s ne vr,count) & test=test+count
ii=where(x_s ne x,count) & test=test+count
ii=where(Tnorm_s ne Tnorm,count) & test=test+count
ii=where(mu_s ne mu,count) & test=test+count
if test le 0 then New_Grid=0
endif
New_Protons=1
if type_of(Ti_s) ne 0 then begin
test=0
ii=where(Ti_s ne Ti,count) & test=test+count
ii=where(n_s ne n,count) & test=test+count
ii=where(vxi_s ne vxi,count) & test=test+count
if test le 0 then New_Protons=0
endif
New_Electrons=1
if type_of(Te_s) ne 0 then begin
test=0
ii=where(Te_s ne Te,count) & test=test+count
ii=where(n_s ne n,count) & test=test+count
if test le 0 then New_Electrons=0
endif
New_fH=1
if type_of(fH_s) ne 0 then begin
ii=where(fH_s ne fH,count)
if count le 0 then New_fH=0
endif
New_Simple_CX=1
if type_of(Simple_CX_s) ne 0 then begin
ii=where(Simple_CX_s ne Simple_CX,count)
if count le 0 then New_Simple_CX=0
endif
New_H2_Seed=1
if type_of(fH2_s) ne 0 then begin
ii=where(fH2_s ne fH2,count)
if count le 0 then New_H2_Seed=0
endif
New_HP_Seed=1
if type_of(nHP_s) ne 0 then begin
test=0
ii=where(nHP_s ne nHP,count) & test=test+count
ii=where(THP_s ne THP,count) & test=test+count
if test le 0 then New_HP_Seed=0
endif
New_ni_correct=1
if type_of(ni_correct_s) ne 0 then begin
ii=where(ni_correct_s ne ni_correct,count)
if count le 0 then New_ni_correct=0
endif
Do_sigv= New_Grid or New_Electrons
Do_fH_moments= (New_Grid or New_fH) and total(fH) gt 0.0
Do_Alpha_CX= (New_Grid or (type_of(Alpha_CX) eq 0) or New_HP_Seed or New_Simple_CX) and H2_HP_CX
;Do_Alpha_CX is updated in fH2_iteration loop
Do_SIG_CX= (New_Grid or (type_of(SIG_CX) eq 0) or New_Simple_CX) and (Simple_CX eq 0) and Do_Alpha_CX
Do_Alpha_H2_H= (New_Grid or (type_of(Alpha_H2_H) eq 0) or New_fH) and H2_H_EL
Do_SIG_H2_H= (New_Grid or (type_of(SIG_H2_H) eq 0)) and Do_Alpha_H2_H
Do_SIG_H2_H2= (New_Grid or (type_of(SIG_H2_H2) eq 0)) and H2_H2_EL
Do_Alpha_H2_P= (New_Grid or (type_of(Alpha_H2_P) eq 0) or New_Protons or New_ni_correct) and H2_P_EL
; Do_Alpha_H2_P is updated in fH2_iteration loop
Do_SIG_H2_P= (New_Grid or (type_of(SIG_H2_P) eq 0)) and Do_Alpha_H2_P
Do_v_v2= (New_Grid or (type_of(v_v2) eq 0)) and (CI_Test or Do_SIG_CX or Do_SIG_H2_H or Do_SIG_H2_H2 or Do_SIG_H2_P)
nH=dblarr(nx)
VxH=dblarr(nx)
TH=dblarr(nx)+1.0
if Do_fH_moments then begin
if debrief gt 1 then print,prompt+'Computing vx and T moments of fH'
;
; Compute x flow velocity and temperature of atomic species
;
for k=0,nx-1 do begin
nH(k)=total(Vr2pidVr*(fH(*,*,k)#dVx))
if nH(k) gt 0 then begin
VxH(k)=vth*total(Vr2pidVr*(fH(*,*,k)#(Vx*dVx)))/nH(k)
for i=0,nvr-1 do vr2vx2_ran2(i,*)=vr(i)^2+(vx-VxH(k)/vth)^2
TH(k)=(mu*mH)*vth2*total(Vr2pidVr*((vr2vx2_ran2*fH(*,*,k))#dVx))/(3*q*nH(k))
endif
endfor
endif
if New_Grid then begin
if debrief gt 1 then print,prompt+'Computing vr2vx2, vr2vx_vxi2, EH2_P'
;
; Magnitude of total normalized v^2 at each mesh point
;
vr2vx2=dblarr(nvr,nvx,nx)
for i=0,nvr-1 do for k=0,nx-1 do vr2vx2(i,*,k)=vr(i)^2+vx^2
;
; Magnitude of total normalized (v-vxi)^2 at each mesh point
;
vr2vx_vxi2=dblarr(nvr,nvx,nx)
for i=0,nvr-1 do for k=0,nx-1 do vr2vx_vxi2(i,*,k)=vr(i)^2+(vx-vxi(k)/vth)^2
;
; Molecular hydrogen ion energy in local rest frame of plasma at each mesh point
;
EH2_P=mH*vr2vx_vxi2*vth2/q
EH2_P=EH2_P > 0.1 &; sigmav_cx does not handle neutral energies below 0.1 eV
EH2_P=EH2_P < 2.0E4 &; sigmav_cx does not handle neutral energies above 20 keV
;
; Compute Maxwellian H2 distribution at the wall temperature
;
fw_Hat=dblarr(nvr,nvx)
;
; NOTE: Molecular ions have 'normalizing temperature' of 2 Tnorm, i.e., in order to
; achieve the same thermal velocity^2, a molecular ion distribution has to have twice the temperature
; as an atomic ion distribution
;
if (total(SH2) gt 0) or (total(PipeDia) gt 0) then begin
if debrief gt 1 then print,prompt+'Computing fw_Hat'
vx_shift=[0.0]
Tmaxwell=[Twall]
mol=2
create_shifted_Maxwellian,vr,vx,Tmaxwell,vx_shift,mu,mol,Tnorm,_Maxwell
fw_hat=_Maxwell(*,*,0)
endif
endif
if New_Protons then begin
;
; Compute Fi_hat
;
if debrief gt 1 then print,prompt+'Computing fi_Hat'
vx_shift=vxi
Tmaxwell=Ti
mol=1
@create_shifted_maxwellian.include
fi_hat=Maxwell
endif
if Do_sigv then begin
if debrief gt 1 then print,prompt+'Computing sigv'
;
; Compute sigmav rates for each reaction and optionally apply
; CR model corrections of Sawada
;
sigv=dblarr(nx,11)
;________________________________________________________________________________
; Reaction R1: e + H2 -> e + H2(+) + e
;________________________________________________________________________________
sigv(*,1)=sigmav_ion_HH(Te)
if Sawada then sigv(*,1)=sigv(*,1)*3.7/2.0
;________________________________________________________________________________
; Reaction R2: e + H2 -> H(1s) + H(1s)
;________________________________________________________________________________
sigv(*,2)=sigmav_H1s_H1s_HH(Te)
if Sawada then begin
;
; Construct table
;
Te_table=alog([5,20,100]) & Ne_Table=alog([1e14,1e17,1e18,1e19,1e20,1e21,1e22])
fctr_Table=fltarr(7,3)
fctr_Table(*,0)=[2.2, 2.2, 2.1, 1.9, 1.2, 1.1, 1.05]/5.3
fctr_Table(*,1)=[5.1, 5.1, 4.3, 3.1, 1.5, 1.25, 1.25]/10.05
fctr_Table(*,2)=[1.3, 1.3, 1.1, 0.8, 0.38, 0.24, 0.22]/2.1
_Te=Te
_Te=_Te > 5
_Te=_Te < 100
_n=n > 1e14
_n=n < 1e22
fctr=Path_Interp_2D(fctr_Table,Ne_Table,Te_table,alog(_n),alog(_Te))
sigv(*,2)=(1.0+fctr)*sigv(*,2)
endif
;________________________________________________________________________________
; Reaction R3: e + H2 -> e + H(1s) + H*(2s)
;________________________________________________________________________________
sigv(*,3)=sigmav_H1s_H2s_HH(Te)
;________________________________________________________________________________
; Reaction R4: e + H2 -> e + p + H(1s)
;________________________________________________________________________________
sigv(*,4)=sigmav_P_H1s_HH(Te)
if Sawada then sigv(*,4)=sigv(*,4)*1.0/0.6
;________________________________________________________________________________
; Reaction R5: e + H2 -> e + H*(2p) + H*(2s)
;________________________________________________________________________________
sigv(*,5)=sigmav_H2p_H2s_HH(Te)
;________________________________________________________________________________
; Reaction R6: e + H2 -> e + H(1s) + H*(n=3)
;________________________________________________________________________________
sigv(*,6)=sigmav_H1s_Hn3_HH(Te)
;________________________________________________________________________________
; Reaction R7: e + H2(+) -> e + p + H(1s)
;________________________________________________________________________________
sigv(*,7)=sigmav_P_H1s_HP(Te)
;________________________________________________________________________________
; Reaction R8: e + H2(+) -> e + p + H*(n=2)
;________________________________________________________________________________
sigv(*,8)=sigmav_p_Hn2_HP(Te)
;________________________________________________________________________________
; Reaction R9: e + H2(+) -> e + p + p + e
;________________________________________________________________________________
sigv(*,9)=sigmav_P_P_HP(Te)
;________________________________________________________________________________
; Reaction R10: e + H2(+) -> e + H(1s) + H*(n>=2)
;________________________________________________________________________________
sigv(*,10)=sigmav_H1s_Hn_HP(Te)
;________________________________________________________________________________
;
; Total H2 destruction rate (normalized by vth) = sum of reactions 1-6
;
alpha_loss=dblarr(nx)
alpha_loss(*)=n*total(sigv(*,1:6),2)/vth
endif
;
;________________________________________________________________________________
; Set up arrays for charge exchange and elastic collision computations, if needed
;________________________________________________________________________________
;
if Do_v_v2 eq 1 then begin
if debrief gt 1 then print,prompt+'Computing v_v2, v_v, vr2_vx2, and vx_vx'
;
; v_v2=(v-v_prime)^2 at each double velocity space mesh point, including theta angle
;
v_v2=dblarr(nvr,nvx,nvr,nvx,ntheta)
;
; vr2_vx2=(vr2 + vr2_prime - 2*vr*vr_prime*cos(theta) - 2*(vx-vx_prime)^2
; at each double velocity space mesh point, including theta angle
;
vr2_vx2=dblarr(nvr,nvx,nvr,nvx,ntheta)
for m=0,ntheta-1 do begin
for l=0,nvx-1 do begin
for k=0,nvr-1 do begin
for i=0,nvr-1 do begin
v_v2(i,*,k,l,m)=Vr(i)^2+Vr(k)^2-2*Vr(i)*Vr(k)*cos_theta(m)+(Vx(*)-Vx(l))^2
vr2_vx2(i,*,k,l,m)=Vr(i)^2+Vr(k)^2-2*Vr(i)*Vr(k)*cos_theta(m)-2*(Vx(*)-Vx(l))^2
endfor
endfor
endfor
endfor
;
; v_v=|v-v_prime| at each double velocity space mesh point, including theta angle
;
v_v=sqrt(v_v2)
;
; vx_vx=(vx-vx_prime) at each double velocity space mesh point
;
vx_vx=dblarr(nvr,nvx,nvr,nvx)
for j=0,nvx-1 do for l=0,nvx-1 do vx_vx(*,j,*,l)=vx(j)-vx(l)
;
; Set Vr'2pidVr'*dVx' for each double velocity space mesh point
;
Vr2pidVrdVx=dblarr(nvr,nvx,nvr,nvx)
for k=0,nvr-1 do Vr2pidVrdVx(*,*,k,*)=Vr2pidVr(k)
for l=0,nvx-1 do Vr2pidVrdVx(*,*,*,l)=Vr2pidVrdVx(*,*,*,l)*dVx(l)
endif
;
if Simple_CX eq 0 and Do_SIG_CX eq 1 then begin
if debrief gt 1 then print,prompt+'Computing SIG_CX'
;________________________________________________________________________________
; Option (A) was selected: Compute SigmaV_CX from sigma directly.
; In preparation, compute SIG_CX for present velocity space grid, if it has not
; already been computed with the present input parameters
;________________________________________________________________________________
;
; Compute sigma_cx * v_v at all possible relative velocities
;
_Sig=dblarr(nvr*nvx*nvr*nvx,ntheta)
_Sig(*)=v_v*sigma_cx_hh(v_v2*(mH*vth2/q))
;
; Set SIG_CX = vr' x Integral{v_v*sigma_cx} over theta=0,2pi times differential velocity space element Vr'2pidVr'*dVx'
;
SIG_CX=dblarr(nvr*nvx,nvr*nvx)
SIG_CX(*)=Vr2pidVrdVx*(_Sig#dtheta)
;
; SIG_CX is now vr' * sigma_cx(v_v) * v_v (intergated over theta) for all possible ([vr,vx],[vr',vx'])
;
endif
;
if Do_SIG_H2_H eq 1 then begin
if debrief gt 1 then print,prompt+'Computing SIG_H2_H'
;________________________________________________________________________________
; Compute SIG_H2_H for present velocity space grid, if it is needed and has not
; already been computed with the present input parameters
;________________________________________________________________________________
;
; Compute sigma_H2_H * v_v at all possible relative velocities
;
_Sig=dblarr(nvr*nvx*nvr*nvx,ntheta)
_Sig(*)=v_v*Sigma_EL_H_HH(v_v2*(0.5*mH*vth2/q))
;
; NOTE: using H energy here for cross-sections tabulated as H->H2
;
; Set SIG_H2_H = vr' x vx_vx x Integral{v_v*sigma_H2_H} over theta=0,2pi times differential velocity space element Vr'2pidVr'*dVx'
;
SIG_H2_H=dblarr(nvr*nvx,nvr*nvx)
SIG_H2_H(*)=Vr2pidVrdVx*vx_vx*(_Sig#dtheta)
;
; SIG_H2_H is now vr' * vx_vx * sigma_H2_H(v_v) * v_v (intergated over theta) for all possible ([vr,vx],[vr',vx'])
;
endif
;
if Do_SIG_H2_P eq 1 then begin
if debrief gt 1 then print,prompt+'Computing SIG_H2_P'
;________________________________________________________________________________
; Compute SIG_H2_P for present velocity space grid, if it is needed and has not
; already been computed with the present input parameters
;________________________________________________________________________________
;
; Compute sigma_H2_P * v_v at all possible relative velocities
;
_Sig=dblarr(nvr*nvx*nvr*nvx,ntheta)
_Sig(*)=v_v*Sigma_EL_P_HH(v_v2*(0.5*mH*vth2/q))
;
; NOTE: using H energy here for cross-sections tabulated as P->H2
;
; Set SIG_H2_P = vr' x vx_vx x Integral{v_v*sigma_H2_P} over theta=0,2pi times differential velocity space element Vr'2pidVr'*dVx'
;
SIG_H2_P=dblarr(nvr*nvx,nvr*nvx)
SIG_H2_P(*)=Vr2pidVrdVx*vx_vx*(_Sig#dtheta)
;
; SIG_H2_P is now vr' * vx_vx * sigma_H2_P(v_v) * v_v (intergated over theta) for all possible ([vr,vx],[vr',vx'])
;
endif
;
if Do_SIG_H2_H2 eq 1 then begin
if debrief gt 1 then print,prompt+'Computing SIG_H2_H2'
;________________________________________________________________________________
; Compute SIG_H2_H2 for present velocity space grid, if it is needed and has not
; already been computed with the present input parameters
;________________________________________________________________________________
;
; Compute sigma_H2_H2 * vr2_vx2 * v_v at all possible relative velocities
;
_Sig=dblarr(nvr*nvx*nvr*nvx,ntheta)
_Sig(*)=vr2_vx2*v_v*Sigma_EL_HH_HH(v_v2*(mH*mu*vth2/q),/VIS)/8.0
;
; Note: For viscosity, the cross section for D -> D is the same function of
; center of mass energy as H -> H.
;
; Set SIG_H2_H2 = vr' x Integral{vr2_vx2*v_v*sigma_H2_H2} over theta=0,2pi times differential velocity space element Vr'2pidVr'*dVx'
;
SIG_H2_H2=dblarr(nvr*nvx,nvr*nvx)
SIG_H2_H2(*)=Vr2pidVrdVx*(_Sig#dtheta)
;
; SIG_H2_H2 is now vr' * sigma_H2_H2(v_v) * vr2_vx2 * v_v (intergated over theta) for all possible ([vr,vx],[vr',vx'])
;
endif
;
;________________________________________________________________________________
; Compute Alpha_H2_H for inputted fH, if it is needed and has not
; already been computed with the present input parameters
;________________________________________________________________________________
;
if Do_Alpha_H2_H eq 1 then begin
if debrief gt 1 then print,prompt+'Computing Alpha_H2_H'
Alpha_H2_H=dblarr(nvr,nvx,nx)
for k=0,nx-1 do begin
Work(*)=fH(*,*,k)
Alpha_H2_H(*,*,k)=SIG_H2_H#Work
endfor
endif
;
;________________________________________________________________________________
; Compute nH2
;________________________________________________________________________________
for k=0,nx-1 do nH2(k)=total(Vr2pidVr*(fH2(*,*,k)#dVx))
;
if New_H2_Seed then begin
MH2_H2_sum=dblarr(nvr,nvx,nx)
Delta_nH2s=1.0
endif
;
; Compute Side-Wall collision rate
;
gamma_wall=dblarr(nvr,nvx,nx)
for k=0,nx-1 do begin
if PipeDia(k) gt 0.0 then begin
for j=0,nvx-1 do gamma_wall(*,j,k)=2*vr/PipeDia(k)
endif
endfor
;
fH2_Iterate:
;
; This is the iteration entry point for fH2, THP and nHP iteration.
; Save 'seed' values for comparison later
;
fH2s=fH2
nH2s=nH2
THPs=THP
nHPs=nHP
;
;________________________________________________________________________________
; Compute Alpha_CX for present THP and nHP, if it is needed and has not
; already been computed with the present parameters
;________________________________________________________________________________
;
if Do_Alpha_CX eq 1 then begin
if debrief gt 1 then print,prompt+'Computing Alpha_CX'
;
; Set Maxwellian Molecular Ion Distribution Function (assumed to be drifting with ion velocity, vxi)
;
vx_shift=vxi
Tmaxwell=THP
mol=2
@create_shifted_maxwellian.include
fHp_hat=Maxwell
;
if Simple_CX then begin
;________________________________________________________________________________
; Option (B): Use maxwellian weighted <sigma v>
;________________________________________________________________________________
;
; THp/mu at each mesh point
;
THp_mu=dblarr(nvr,nvx,nx)
for k=0,nx-1 do THp_mu(*,*,k)=THp(k)/mu
;
; Molecular Charge Exchange sink rate
;
alpha_cx=sigmav_cx_HH(THp_mu,EH2_P)/vth
for k=0,nx-1 do alpha_cx(*,*,k)=alpha_cx(*,*,k)*nHp(k)
;________________________________________________________________________________
;
endif else begin
;________________________________________________________________________________
; Option (A): Compute SigmaV_CX from sigma directly via SIG_CX
;________________________________________________________________________________
;
alpha_cx=dblarr(nvr,nvx,nx)
for k=0,nx-1 do begin
Work(*)=fHp_hat(*,*,k)*nHp(k)
alpha_cx(*,*,k)=SIG_CX#Work
endfor
if do_alpha_cx_test then begin
alpha_cx_test=sigmav_cx_HH(THp_mu,EH2_P)/vth
for k=0,nx-1 do alpha_cx_test(*,*,k)=alpha_cx_test(*,*,k)*nHp(k)
print,'Compare alpha_cx and alpha_cx_test'
press_return
endif
endelse
endif
;________________________________________________________________________________
; Compute Alpha_H2_P for present Ti and ni (optionally correcting for nHP),
; if it is needed and has not already been computed with the present parameters
;________________________________________________________________________________
if Do_Alpha_H2_P eq 1 then begin
if debrief gt 1 then print,prompt+'Computing Alpha_H2_P'
Alpha_H2_P=dblarr(nvr,nvx,nx)
ni=n
if ni_correct then ni=(n-nHp) > 0.0
for k=0,nx-1 do begin