0001 function equil = imake_equil_jd(dkepath,basestr,signs,external,workdir,opt_gui,select)
0002
0003
0004
0005 if nargin < 7,
0006 select = struct;
0007 end
0008
0009 if nargin < 6,
0010 opt_gui = true;
0011 end
0012
0013 if nargin < 5,
0014 workdir = pwd;
0015 end
0016
0017 if nargin < 4,
0018 external = '';
0019 end
0020
0021 if nargin < 3,
0022 signs = '';
0023 end
0024
0025 if ~isfield(select,'style'),
0026 select.style = struct;
0027 end
0028
0029 select.style.textwidth = 420;
0030
0031 eqdsk_path = './';
0032
0033 equil = '';
0034 Zaxis = 0;
0035 Zeff = NaN;
0036
0037 e_opt = 0;
0038 f_opt = 0;
0039 type = NaN;
0040
0041
0042
0043 disp(' ')
0044 disp('*********** MAGNETIC EQUILIBRIUM ***********')
0045 disp(' ')
0046
0047 if isfield(select,'i_opt'),
0048 i_opt = select.i_opt;
0049 elseif isfield(external,'magnetic'),
0050
0051 titletext = 'For the MAGNETIC equilibrium :';
0052 s = {'Use data loaded from the external database',...
0053 'Use data from an existing equilibrium file',...
0054 'Build the magnetic equilibrium data locally'};
0055
0056 i_opt = iselect_jd(s,titletext,opt_gui,select.style,1) - 1;
0057
0058 else
0059
0060 titletext = 'For the MAGNETIC equilibrium :';
0061 s = {'Use data from an existing equilibrium file',...
0062 'Build the magnetic equilibrium data locally'};
0063
0064 i_opt = iselect_jd(s,titletext,opt_gui,select.style,1);
0065
0066 end
0067
0068 if i_opt == 0,
0069
0070 shotnum = external.shotnum;
0071 shotime = external.shotime;
0072 if ischar(shotime),
0073 shotimestr = shotime;
0074 shotime = str2double(shotimestr);
0075 else
0076 shotimestr = num2str(shotime,'%6.4f');
0077 end
0078
0079 if strcmp(basestr,'TCV') || strcmp(basestr,'EAST') || strcmp(basestr,'HL2A') || strcmp(basestr,'CMOD'),
0080
0081 if isfield(external.magnetic,'eqdsk'),
0082
0083 eqdsk = external.magnetic.eqdsk;
0084
0085 elseif isfield(external.magnetic,'eqdsk_path') && isfield(external.magnetic,'eqdsk_file')
0086
0087 eqdsk_path = external.magnetic.eqdsk_path;
0088 eqdsk_file = external.magnetic.eqdsk_file;
0089
0090 if ~exist(eqdsk_path,'dir'),
0091 disp(['WARNING : The directory ',eqdsk_path,' does not exist. local directory enforced.'])
0092 eqdsk_path = './';
0093 end
0094 if ~exist([eqdsk_path,eqdsk_file],'file'),
0095 error(['The file ',eqdsk_path,eqdsk_file,' does not exist. Aborted.'])
0096 end
0097
0098 eqdsk = fileread([eqdsk_path,eqdsk_file]);
0099
0100 end
0101
0102 e_opt = -1;
0103
0104 elseif strcmp(basestr,'TS'),
0105
0106 equil_magnetic_data = external.magnetic;
0107
0108 elseif strcmp(basestr,'COMPASS'),
0109
0110 equil_magnetic_data = external.magnetic;
0111
0112 eqdsk.nnx = length(external.magnetic.efit.R);
0113 eqdsk.nny = length(external.magnetic.efit.Z);
0114 eqdsk.nnr = length(external.magnetic.efit.psi_n);
0115
0116 eqdsk.xbox = external.magnetic.efit.R(end) - external.magnetic.efit.R(1);
0117 eqdsk.ybox = abs(external.magnetic.efit.Z(end) - external.magnetic.efit.Z(1));
0118 eqdsk.Rv = external.magnetic.efit.R_geom_axis;
0119 eqdsk.Rmin = external.magnetic.efit.R(1);
0120 eqdsk.Zv = 0;
0121
0122 eqdsk.Rp = external.magnetic.efit.R_mag_axis;
0123 eqdsk.Zp = external.magnetic.efit.Z_mag_axis;
0124 eqdsk.psip = external.magnetic.efit.psi_mag_axis;
0125 eqdsk.psia = external.magnetic.efit.psi_lcfs;
0126 eqdsk.Bv = external.magnetic.efit.B_vac_R_geom*external.magnetic.efit.R_geom_axis/external.magnetic.efit.R_mag_axis;
0127
0128 eqdsk.Ip = external.magnetic.efit.Ip;
0129
0130 eqdsk.rF = external.magnetic.efit.RBphi.';
0131 eqdsk.rp = external.magnetic.efit.p_stat;
0132 eqdsk.rFdFdpsi = external.magnetic.efit.FFprime;
0133 eqdsk.rdpdpsi = external.magnetic.efit.pprime_static;
0134
0135 eqdsk.xypsi = external.magnetic.efit.psi_RZ;
0136
0137 eqdsk.rq = external.magnetic.efit.q.';
0138
0139 e_opt = -1;
0140
0141 elseif strcmp(basestr,'JET'),
0142
0143 equil_magnetic_data = external.magnetic;
0144
0145 e_opt = -3;
0146
0147 elseif strcmp(basestr,'AUG'),
0148
0149 equil_magnetic_data = external.magnetic;
0150
0151 e_opt = -4;
0152
0153 else
0154
0155 disp('-----> function not yet implemented for this machine.')
0156 return
0157 end
0158
0159 elseif i_opt == 1,
0160
0161 d_opt = input_dke_yp('Format of the magnetic equilibrium data, 0 (eqdsk), 1 (itm), (2) (AMR): ',0,{0,1,2},'Invalid option');
0162
0163 if d_opt == 0,
0164
0165
0166
0167 [eqdsk_file,eqdsk_path] = igetfile_jd(opt_gui,{'eqdsk*.*';'EQDSK*.*';'*.*'},'Please select the EQDSK file',workdir);
0168
0169 if eqdsk_file == 0,
0170 disp('-----> No eqdsk file selected.')
0171 return
0172 end
0173
0174 eqdsk = fileread([eqdsk_path,eqdsk_file]);
0175
0176 it = strfind(eqdsk_file,'t');
0177
0178 if (~isempty(strfind(eqdsk_file,'EQDSK.')) || ~isempty(strfind(eqdsk_file,'eqdsk.'))) && ~isempty(it),
0179 shotnum = eqdsk_file(7:it-1);
0180 shotimestr = eqdsk_file(it+1:end);
0181 else
0182 shotnum = '0';
0183 shotimestr = '0';
0184 end
0185
0186 shotnum_new = input(['Please provide the shot number [',shotnum,'] : '],'s');
0187 if ~isempty(shotnum_new),
0188 shotnum = shotnum_new;
0189 elseif isempty(shotnum),
0190 shotnum = '0';
0191 end
0192
0193 shotime_new = input(['Please provide the shot time [',shotimestr,'] : '],'s');
0194 if ~isempty(shotime_new),
0195 shotimestr = shotime_new;
0196 elseif isempty(shotimestr),
0197 shotimestr = '0';
0198 end
0199
0200 shotime = str2double(shotimestr);
0201
0202 e_opt = -1;
0203
0204 elseif d_opt == 1,
0205
0206 disp('-----> loading of ITM equilibrium not yet implemented.')
0207 return
0208
0209 elseif d_opt == 2,
0210
0211
0212
0213 [amr_psi_file,amr_psi_path] = igetfile_jd(opt_gui,{'m_profs_*_Psi.dat';'*.*'},'Please select the AMR equilibrium file containing the Psi vector',workdir);
0214
0215 if amr_psi_file == 0,
0216 disp('-----> No AMR equilibrium file selected.')
0217 return
0218 end
0219
0220 if (strcmp(amr_psi_file(1:8),'m_profs_')),
0221 it1 = 8 + strfind(amr_psi_file(9:end),'_');
0222 shotnum = amr_psi_file(9:it1(1)-1);
0223 it2 = it1(1) + strfind(amr_psi_file(it1+1:end),'_');
0224 shotimestr = amr_psi_file(it1(1)+1:it2(1)-1);
0225 else
0226 disp('-----> The AMR equilibrium file selected is not valid.')
0227 return
0228 end
0229
0230 e_opt = -2;
0231 n_opt = 0;
0232
0233 shotime = str2double(shotimestr);
0234
0235 end
0236
0237 end
0238
0239 if e_opt == 0,
0240
0241 e_opt = input_dke_yp('Select magnetic equilibrium solver, 0 (analytic circular concentric), 1 (helena) : ',0,{0,1},'Invalid option');
0242
0243 end
0244
0245 if i_opt == 2,
0246
0247 [equil_magnetic_data,shotnum,shotimestr] = imake_equil_ideal(e_opt);
0248
0249 shotime = str2double(shotimestr);
0250
0251 end
0252
0253
0254
0255 disp(' ')
0256 disp('*********** EQUILIBRIUM PROFILES ***********')
0257 disp(' ')
0258
0259 if e_opt ~= -2,
0260 if isfield(select,'i_opt'),
0261 i_opt = select.i_opt;
0262 elseif isfield(external,'prof'),
0263
0264 titletext = 'For the equilibrium PROFILES:';
0265 s = {'Use data loaded from the external database',...
0266 'Use data from an existing equilibrium profiles file',...
0267 'Build the equilibrium profiles data locally'};
0268
0269 i_opt = iselect_jd(s,titletext,opt_gui,select.style,1) - 1;
0270
0271 else
0272
0273 titletext = 'For the equilibrium PROFILES:';
0274 s = {'Use data from an existing equilibrium profiles file',...
0275 'Build the equilibrium profiles data locally'};
0276
0277 i_opt = iselect_jd(s,titletext,opt_gui,select.style,1);
0278
0279 end
0280
0281 if i_opt == 0,
0282
0283 equil_prof_data = external.prof;
0284
0285 if strcmp(basestr,'TCV'),
0286
0287 if isfield(equil_prof_data,'shot'),
0288 type = 1;
0289 else
0290 type = -2;
0291
0292 datamod.Zaxis = equil_prof_data.zaxis;
0293 Ip = - equil_prof_data.ip;
0294 signs.Bt = - equil_prof_data.btsign;
0295 signs.Ip = - equil_prof_data.ipsign;
0296
0297 profname = rmfield(equil_prof_data,{'zaxis','ip','btsign','ipsign'});
0298
0299 end
0300
0301 n_opt = 0;
0302
0303 elseif strcmp(basestr,'TS'),
0304
0305 n_opt = -1;
0306
0307 elseif strcmp(basestr,'JET'),
0308
0309 if isempty(equil_prof_data.pne) || isempty(equil_prof_data.pTe),
0310
0311 f_opt = 1;
0312 n_opt = -1;
0313
0314 else
0315
0316 type = 1.2;
0317 n_opt = 0;
0318
0319 end
0320 elseif strcmp(basestr,'EAST') || strcmp(basestr,'HL2A'),
0321
0322 type = 1.4;
0323 n_opt = 0;
0324
0325 elseif strcmp(basestr,'COMPASS'),
0326
0327 type = 1.6;
0328 n_opt = 0;
0329
0330 elseif strcmp(basestr,'AUG'),
0331
0332 type = -2;
0333 profname = equil_prof_data;
0334 n_opt = 0;
0335
0336 if isfield(equil_prof_data,'Ip')
0337 Ip = equil_prof_data.Ip;
0338 end
0339
0340 signs = struct;
0341 else
0342
0343 disp('-----> function not yet implemented for this machine.')
0344 return
0345 end
0346
0347 elseif i_opt == 1,
0348
0349 n_opt = input_dke_yp(['Option for equilibrium profiles file : \n',...
0350 '0 (existing data file)\n',...
0351 '1 (ART profiles)\n',...
0352 '2 (TORAY profiles)\n',...
0353 '3 (EVE profiles)\n'],...
0354 0,{0,1,2,3},'Invalid selection');
0355
0356 [data_file,data_path] = igetfile_jd(opt_gui,'*.*',['Please select the data file for shot ',shotnum,' at time ',shotimestr,' :'],eqdsk_path);
0357
0358 if data_file == 0,
0359 return
0360 end
0361
0362 profname = [data_path,data_file];
0363
0364 if n_opt == 0,
0365
0366 if strcmp(basestr,'TCV'),
0367
0368 type = 1;
0369
0370 equil_prof_data = load([data_path,data_file],'shot','rho','ne','te','ti','zeff','zaxis','btsign','ipsign','Ip');
0371
0372 elseif strcmp(basestr,'TS'),
0373
0374 type = -1;
0375
0376 elseif strcmp(basestr,'ITER'),
0377
0378 type = 0.2;
0379
0380 elseif strcmp(basestr,'EAST') || strcmp(basestr,'HL2A'),
0381
0382 type = 1.4;
0383
0384 elseif strcmp(basestr,'COMPASS'),
0385
0386 type = 1.6;
0387
0388 else
0389 disp(['-----> profiles data files not treated yet for :',basestr]);
0390 return
0391 end
0392
0393 elseif n_opt == 1,
0394
0395 type = 1;
0396
0397 load([data_path,data_file],'ARTsimprof');
0398
0399 equil_prof_data.rho = ARTsimprof.rho.data.';
0400 equil_prof_data.ne = ARTsimprof.n.data.'*1e20;
0401 equil_prof_data.te = ARTsimprof.X.data.';
0402 equil_prof_data.ti = NaN;
0403
0404 equil_prof_data.zeff = NaN;
0405 equil_prof_data.zaxis = NaN;
0406
0407 equil_prof_data.shot = str2num(data_file(4:8));
0408
0409 elseif n_opt == 2,
0410
0411 disp('-----> profiles data files not treated yet for TORAY output');
0412 return
0413
0414 elseif n_opt == 3,
0415
0416 type = 1.3;
0417
0418 end
0419
0420 elseif i_opt == 2,
0421
0422 equil_prof_data = imake_equil_ideal;
0423
0424 n_opt = -1;
0425
0426 end
0427 end
0428
0429
0430
0431 shotid = [basestr,'_',shotnum,'_',shotimestr];
0432
0433 if ~isfield(select,'idequil') || ~select.idequil
0434 id_equil = input(['Please provide a name for the equilibrium file [',shotid,'] : '],'s');
0435 else
0436 id_equil = '';
0437 end
0438
0439 if isempty(id_equil),
0440 id_equil = shotid;
0441 end
0442
0443 if (~isfield(select,'disp') || ~select.disp) && strcmp(input_dke_yp('Do you want to display the equilibrium','n',{'y','n'},'',[1,1]),'y'),
0444 display_mode = 2;
0445 else
0446 display_mode = 1;
0447 end
0448
0449 if (~isfield(select,'print') || ~select.print)
0450 p_opt = input_dke_yp('Option for equilibrium figures? -1 (nothing), 0 (print), 1 (print and save), 2 (save)',-1,-1:2,'',[1,1]);
0451 else
0452 p_opt = -1;
0453 end
0454
0455 dispname = '';
0456 savename = id_equil;
0457
0458 if type == 1 && n_opt <= 1,
0459
0460 ishot = find(equil_prof_data.shot == str2num(shotnum));
0461
0462 if isempty(ishot),
0463 disp(['-----> Profiles data not found for shot :',shotnum]);
0464 return
0465 end
0466
0467 if length(ishot) > 1,
0468 disp(['-----> WARNING : Several entries are found for shot :',shotnum,', first entry selected']);
0469 ishot = ishot(1);
0470 end
0471
0472 prho = equil_prof_data.rho(ishot,:);
0473 pne = equil_prof_data.ne(ishot,:);
0474 pTe = equil_prof_data.te(ishot,:);
0475
0476 datamod.Zaxis = equil_prof_data.zaxis(ishot);
0477
0478 if any(isnan(equil_prof_data.ti(ishot,:))),
0479 datamod.TeTiratio = 6;
0480 end
0481
0482 if isscalar(equil_prof_data.zeff(ishot,:)),
0483 if ~isnan(equil_prof_data.zeff(ishot,:)),
0484 datamod.Zeff = equil_prof_data.zeff(ishot,:);
0485 else
0486 datamod.Zeff = 3;
0487 end
0488 end
0489
0490 if isfield(equil_prof_data,'Ip')
0491 datamod.Ip = - equil_prof_data.Ip(ishot)/1e6;
0492 else
0493 datamod.Ip = NaN;
0494 end
0495
0496 s = struct;
0497
0498 s.expert.Zaxis = 0;
0499 s.expert.TeTiratio = 0;
0500 s.expert.Zeff = 0;
0501 s.expert.Ip = 0;
0502
0503 datamod = imod_struct_jd(datamod,'equil.data',1,opt_gui,s,select.style);
0504
0505 if isfield(datamod,'TeTiratio')
0506 pTi = pTe/datamod.TeTiratio;
0507 else
0508 pTi = equil_prof_data.ti(ishot,:);
0509 end
0510
0511 if isfield(datamod,'Zeff'),
0512 pZeff = datamod.Zeff*ones(size(prho));
0513 else
0514 pZeff = equil_prof_data.zeff(ishot,:);
0515 end
0516
0517 Ip = datamod.Ip*1e6;
0518
0519 profname = [id_equil,'.prof'];
0520
0521
0522
0523 if ~isstruct(signs)
0524 if exist('btsign','var') && exist('ipsign','var'),
0525 signs.Bt = - btsign;
0526 signs.Ip = - ipsign;
0527 end
0528
0529 if isfield(equil_prof_data,'btsign') && isfield(equil_prof_data,'ipsign'),
0530 signs.Bt = - equil_prof_data.btsign;
0531 signs.Ip = - equil_prof_data.ipsign;
0532 end
0533 end
0534
0535 save(profname,'prho','pne','pTe','pTi','pZeff','-mat');
0536
0537 elseif type == 1.2 && n_opt == 0,
0538
0539 Zeff = equil_prof_data.Zeff;
0540 pspsin = sqrt(equil_prof_data.psin);
0541 pne = equil_prof_data.pne;
0542 pTe = equil_prof_data.pTe/1e3;
0543
0544 TeTi = input_dke_yp('No data on Ti was found. Please provide an estimated value for Te/Ti',1,'','',[1,1]);
0545
0546 pTi = pTe/TeTi;
0547
0548 profname = [id_equil,'.prof'];
0549
0550 save(profname,'pspsin','pne','pTe','pTi','Zeff','-mat');
0551
0552 elseif type == 1.4 && n_opt == 0,
0553
0554 profname = [id_equil,'.prof'];
0555
0556 Zeff = equil_prof_data.species.Zeff;
0557 TeTi = equil_prof_data.TeTi;
0558 frac_H = equil_prof_data.species.frac_H;
0559 frac_D = equil_prof_data.species.frac_D;
0560 Zimp = equil_prof_data.species.Zimp;
0561 Aimp = equil_prof_data.species.Aimp;
0562
0563 if isfield(external,'thomson') == 1 & isfield(equil_prof_data,'rho') == 1,
0564 thomson_mode = input_dke_yp('Do you want to calculate Te and ne profiles from raw Thomson scattering measurements','n',{'y','n'},'',[1,1]);
0565 elseif isfield(external,'thomson') == 1 & isfield(equil_prof_data,'rho') == 0,
0566 thomson_mode = 1;
0567 elseif isfield(external,'thomson') == 0 & isfield(equil_prof_data,'rho') == 1,
0568 thomson_mode = 0;
0569 end
0570
0571 if thomson_mode == 1
0572 save(profname,'frac_H','frac_D','Zimp','Aimp','Zeff','-mat');
0573 else
0574 prho = equil_prof_data.rho;
0575 pne = equil_prof_data.ne;
0576 pTe = equil_prof_data.Te;
0577 pTi = pTe/TeTi;
0578
0579 save(profname,'prho','pne','pTe','pTi','frac_H','frac_D','Zimp','Aimp','Zeff','-mat');
0580
0581 end
0582
0583
0584
0585 if ~isstruct(signs)
0586 if isfield(external,'jbdirections'),
0587 if strcmp(external.jbdirections.Ip,'clockwise')
0588 signs.Ip = 1;
0589 else
0590 signs.Ip = -1;
0591 end
0592
0593 if strcmp(external.jbdirections.Bt,'clockwise')
0594 signs.Bt = 1;
0595 else
0596 signs.Bt = -1;
0597 end
0598 end
0599 end
0600
0601 elseif type == 1.6 && n_opt == 0,
0602
0603 profname = [id_equil,'.prof'];
0604
0605 pspsin = equil_prof_data.pspsin;
0606 pne = equil_prof_data.pne.';
0607 pTe = equil_prof_data.pTe.';
0608 pTi = equil_prof_data.pTi.';
0609 Zeff = equil_prof_data.Zeff;
0610 zZi = equil_prof_data.zZi;
0611 zmi = equil_prof_data.zmi;
0612 fi = equil_prof_data.fi;
0613
0614 save(profname,'pspsin','pne','pTe','pTi','Zeff','zZi','zmi','fi','-mat');
0615
0616 end
0617
0618 if e_opt >= 0 || (e_opt <= -2 && e_opt >= -3),
0619
0620 np = input_dke_yp('Number of equilibrium radial points',101,1:2:1001,'np must be an odd number',[1,1]);
0621
0622 nt = input_dke_yp('Number of equilibrium poloidal points',65,3:2:1001,'nt must be an odd number',[1,1]);
0623
0624 end
0625
0626 if e_opt >= 0,
0627
0628
0629
0630 equil_magnetic.Rp = equil_magnetic_data.Rp;
0631 equil_magnetic.Zp = equil_magnetic_data.Zp;
0632
0633 [mpsin,equil_magnetic.psi_apRp,equil_magnetic.theta,equil_magnetic.ptx,equil_magnetic.pty,equil_magnetic.ptBx,equil_magnetic.ptBy,equil_magnetic.ptBPHI,mBpp,mq_Rpap,mj,mmag,Ip_test,mrho,mspsin,mspsinT] = ...
0634 idealequilcyl_yp(equil_magnetic_data.ap,equil_magnetic_data.Rp,equil_magnetic_data.Zp,equil_magnetic_data.Bt,equil_magnetic_data.Ip,equil_magnetic_data.qmin_Rpap,equil_magnetic_data.eq,equil_magnetic_data.qopt,np,nt);
0635
0636 elseif e_opt == -1,
0637
0638 if isfield(select,'busyhandle'),
0639 set(select.busyhandle,'Visible','on');drawnow;
0640 end
0641 [equil_magnetic,mrho,mspsin,mspsinT,mp_magnetic] = equil_magnetic_EFIT_jd(eqdsk,display_mode,p_opt,dispname,savename);
0642
0643 if isfield(select,'busyhandle'),
0644 set(select.busyhandle,'Visible','off');drawnow;
0645 end
0646
0647 [np,nt] = size(equil_magnetic.ptx);
0648
0649
0650
0651 if ~exist('datamod','var'),
0652 Zp_offset = iselect_jd({},'Vertical equilibrium offset (vacuum)',opt_gui,select.style,Zaxis,'edit');
0653 if isempty(Zp_offset) || ~isnumeric(Zp_offset),
0654 Zp_offset = Zaxis;
0655 end
0656 equil_magnetic.Zp = equil_magnetic.Zp + Zp_offset;
0657
0658 else
0659 equil_magnetic.Zp = datamod.Zaxis;
0660 end
0661
0662 elseif e_opt == -2,
0663
0664 [equil_magnetic,equil_prof] = make_equil_AMR_jd(amr_psi_path,shotnum,shotimestr,display_mode,p_opt,dispname,savename,np,nt);
0665
0666 elseif e_opt == -3,
0667
0668 [equil_magnetic,mrho,mspsin,mspsinT] = equil_magnetic_rz2pt_jd(equil_magnetic_data,display_mode,p_opt,dispname,savename,np,nt);
0669
0670 elseif e_opt == -4,
0671
0672 [equil_magnetic,mrho,mspsin] = equil_magnetic_AUG_jd(equil_magnetic_data);
0673 mspsinT = [];
0674
0675 end
0676
0677 if f_opt == 1,
0678
0679 equil_prof_data = fitprof_ideal_jd(equil_prof_data);
0680
0681 end
0682
0683 if n_opt == -1,
0684
0685 [equil_prof_data.prho,equil_prof.pTe,equil_prof.pne,equil_prof.pzTi,equil_prof.pzni,equil_prof.zZi,equil_prof.zmi,equil_prof_data.fi,equil_prof_data.pkin] = ...
0686 idealprof_yp(equil_prof_data.Zi,equil_prof_data.mi,equil_prof_data.fi,equil_prof_data.Te0,equil_prof_data.Tea,equil_prof_data.eTe1,equil_prof_data.eTe2,equil_prof_data.ne0,equil_prof_data.nea,equil_prof_data.ene1,equil_prof_data.ene2,equil_prof_data.Ti0,equil_prof_data.Tia,equil_prof_data.eTi1,equil_prof_data.eTi2,equil_prof_data.Zeff0,equil_prof_data.Zeffa,equil_prof_data.eZeff,equil_prof_data.xZeff,np);
0687
0688 end
0689
0690 if e_opt == 1,
0691
0692
0693
0694 mj_min = min(mj);
0695 if mj_min < 0,
0696 mj = mj - mj_min;
0697 end
0698
0699 [equil_magnetic,mrho,mspsin,mspsinT,mp_magnetic] = equil_magnetic_helmex77_yp(equil_magnetic_data.ap,equil_magnetic_data.Rp,equil_magnetic_data.Zp,equil_magnetic_data.xpoint,equil_magnetic_data.Bt,equil_magnetic_data.Ip,equil_magnetic.psi_apRp*equil_magnetic_data.Rp/equil_magnetic_data.ap,mj,equil_prof_data.pkin,nt,display_mode);
0700
0701 end
0702
0703 if n_opt >= 0 && e_opt ~= -2,
0704
0705 if strcmp(basestr,'EAST') & (thomson_mode == 1 | thomson_mode == 'y')
0706
0707 R_hcn = 1.9;
0708 [equil_magnetic_nete,ne_lav] = EAST_thomson_scattering_Tene_fit(equil_magnetic,external.thomson.prof.R_ts,external.thomson.prof.Z_ts,external.thomson.prof.Te_ts,external.thomson.prof.dTe_ts,external.thomson.prof.ne_ts,external.thomson.prof.dne_ts,R_hcn,display_mode-1);
0709
0710 disp('--> Te and ne profiles determined from raw Thomson scattering measurements.');
0711
0712 prho = mrho;
0713 pTe = equil_magnetic_nete.pTe';
0714 pne = equil_magnetic_nete.pne';
0715
0716 pTi = pTe/TeTi;
0717
0718 save([savename,'.prof'],'prho','pne','pTe','pTi','ne_lav','-append');
0719
0720 end
0721
0722 equil_prof = equil_prof_jd(profname,type,mrho,mspsin,mspsinT,display_mode,p_opt,dispname,savename,Zeff);
0723
0724 if exist('ne_lav')
0725 equil_prof.ne_lav = ne_lav;
0726 else
0727 equil_prof.ne_lav = '';
0728 end
0729
0730 end
0731
0732 if exist('mp_magnetic','var'),
0733
0734 if (~isfield(select,'comp') || ~select.comp) && strcmp(input_dke_yp('Do you want to compare the pressure profiles','n',{'y','n'},'',[1,1]),'y'),
0735
0736
0737
0738 qe = pc_dke_yp;
0739
0740 mp_profiles = (equil_prof.pTe.*equil_prof.pne + sum(equil_prof.pzTi.*equil_prof.pzni,1))*qe*1e3;
0741
0742 figure(14),clf
0743
0744 graph1D_jd(mrho,mp_profiles/1e5,0,0,'r/a','p (Bar)','',NaN,[0,1],NaN,'-','none','b',2,20,gca,0.9,0.7,0.7);
0745 graph1D_jd(mrho,mp_magnetic/1e5,0,0,'','','',{'T_en_e + T_in_i','Equilibrium code'},[0,1],NaN,'-','none','r',2);
0746 set(gca,'XTick',0:0.2:1);
0747 print_jd(p_opt,[savename,'_p_comp']);
0748 end
0749
0750 end
0751
0752
0753
0754 if ~isstruct(signs),
0755 s_opt = input_dke_yp('Do you want to impose the signs of Bt and Ip [y/n]','n',{'y','n'});
0756
0757 if strcmp(s_opt,'y'),
0758 signs.Bt = input_dke_yp('Signs of Bt (positive if clockwise from top)',sign(equil_magnetic.ptBPHI(1,1)),{-1,1});
0759 signs.Ip = input_dke_yp('Signs of Ip (positive if clockwise from top)',sign(equil_magnetic.psi_apRp(end)),{-1,1});
0760 else
0761 signs = struct;
0762 end
0763 end
0764 if ~isfield(signs,'Ip'),
0765 signs.Ip = sign(equil_magnetic.psi_apRp(end));
0766 end
0767 if ~isfield(signs,'Bt'),
0768 signs.Bt = sign(equil_magnetic.ptBPHI(1,1));
0769 end
0770
0771 if signs.Ip*equil_magnetic.psi_apRp(end) < 0,
0772 equil_magnetic.psi_apRp = -equil_magnetic.psi_apRp;
0773 equil_magnetic.ptBx = -equil_magnetic.ptBx;
0774 equil_magnetic.ptBy = -equil_magnetic.ptBy;
0775 end
0776
0777 if signs.Bt*equil_magnetic.ptBPHI(1,1) < 0,
0778 equil_magnetic.ptBPHI = -equil_magnetic.ptBPHI;
0779 end
0780
0781
0782
0783 equil = conc_struct_jd(equil_magnetic,equil_prof);
0784
0785 equil.signs = signs;
0786 equil.id = id_equil;
0787 equil.tokamak = basestr;
0788 equil.shotnum = shotnum;
0789 equil.shotime = shotime;
0790
0791 if exist('Ip','var'),
0792 equil.ip = Ip;
0793 else
0794 equil.ip = NaN;
0795 end
0796
0797 if e_opt == -1 && exist('eqdsk_file','var'),
0798 equil.info.eqdsk = [eqdsk_path,eqdsk_file];
0799 elseif e_opt == -2,
0800 equil.info.datapath = amr_psi_path;
0801 end
0802
0803
0804 if (~isfield(select,'consist') || ~select.consist) && strcmp(input_dke_yp('Do you want to check the equilibrium consistency','n',{'y','n'},'',[1,1]),'y'),
0805 equil = equilconsistency_yp(equil,'',1);
0806 end
0807
0808 equil.info.txt = {...
0809 ['Plasma major radius (m) : ',num2str(equil.Rp)],...
0810 ['Plasma vertical position (m) : ',num2str(equil.Zp)],...
0811 ['Plasma minor radius (m) : ',num2str(equil.ptx(end,1))],...
0812 ['Magnetic field on axis (T) : ',num2str(equil.ptBPHI(1,1))],...
0813 ['Poloidal flux on LCFS (W) : ',num2str(equil.psi_apRp(end)*2*pi*equil.Rp/equil.ptx(end,1))],...
0814 ['Plasma current (MA) : ',num2str(equil.ip/1e6)],...
0815 ['Central electron temperature (keV) : ',num2str(equil.pTe(1))],...
0816 ['Central electron density (10^19 m^-3) : ',num2str(equil.pne(1)/1e19)]};
0817
0818
0819
0820