Teneprof_jd

 This function reads the temperature and density parameters from text or
 binary files in which the columns are respectively the radial position r/a, Te,
 ne, Ti and ni

 by J. Decker (jodecke@alum.mit.edu)  and Y. Peysson (yves.peysson@cea.fr)

0001 function [prho,pspsin,pspsinT,pTe,pne,pzTi,pzni,zZi,zmi] = Teneprof_jd(filename,type,display_mode,p_opt,Zeff,dispname,savename)
0002 %
0003 % This function reads the temperature and density parameters from text or
0004 % binary files in which the columns are respectively the radial position r/a, Te,
0005 % ne, Ti and ni
0006 %
0007 % by J. Decker (jodecke@alum.mit.edu)  and Y. Peysson (yves.peysson@cea.fr)
0008 %
0009 if nargin < 6
0010     savename = filename;
0011 end
0012 if nargin < 5
0013     dispname = filename;
0014 end
0015 %
0016 prho = [];
0017 pspsin = [];
0018 pspsinT = [];
0019 %
0020 if type == -2,% filename is a structure - TCV,AUG
0021     pTe = filename.te.data.';
0022     prho.pTe = filename.te.rho.';
0023     pne = filename.ne.data.';
0024     prho.pne = filename.ne.rho.';
0025     if isfield(filename.ti,'data'),
0026         pTi = filename.ti.data.';
0027         prho.pzTi = filename.ti.rho.';
0028     else
0029         if isscalar(filename.ti),%AUG
0030             pTi = pTe*filename.ti;
0031         else
0032             pTi = pTe/6;%%%% default for TCV
0033         end
0034         prho.pzTi = prho.pTe;
0035     end
0036     if isfield(filename.zeff,'data'),
0037         pZeff = interp1(filename.zeff.rho.',filename.zeff.data.',prho.pne,'linear','extrap');%interpolate on ne grid for zni calculation
0038     elseif ~isnan(filename.zeff),
0039         pZeff = filename.zeff;
0040     else
0041         pZeff = 3;%%%% default for TCV
0042     end
0043     prho.pzni = prho.pne;
0044 elseif type == -1,
0045     eval(['load -mat ',filename,' prho pne pTe pzTi pzni zZi zmi']);
0046 elseif type == 0,%ascii
0047     [prho,pTe,pne,pTi,pni] = textread(filename,'%n%n%n%n%n');
0048 elseif type == 0.1,
0049     [psi,pTe,pne,pZeff] = textread(filename,'%n%n%n%n');    
0050     pTi = pTe;
0051 elseif type == 0.2,
0052     [pspsinT,pTe,pTi,pne,pnH,pnalpha,pnimp,pZeff,dummy] = textread(filename,'%n%n%n%n%n%n%n%n%n');    
0053 elseif type == 1,%.mat
0054     eval(['load -mat ',filename,' prho pTe pne pTi pZeff']);
0055 elseif type == 1.1 ,%.mat
0056     eval(['load -mat ',filename,' prho pTe pne pTi Zeff']);
0057 elseif type == 1.2 ,%.mat
0058     eval(['load -mat ',filename,' pspsin pTe pne pTi Zeff']);    
0059 elseif type == 1.3 ,%.mat, EVE profiles
0060     %
0061     load(filename,'eve','-mat');
0062     %
0063     pspsinT = eve.sg.';
0064     pTe = eve.Ts(:,1).';
0065     pne = eve.ns(:,1).';
0066     pzTi = eve.Ts(:,2:end).';
0067     pzni = eve.ns(:,2:end).';
0068     zZi = eve.Zs(2:end).';
0069     zmi = eve.As(2:end).';
0070     %
0071 elseif type == 1.4 ,%.mat, EAST profiles
0072     %
0073     eval(['load -mat ',filename,' prho pne pTe pTi frac_H frac_D Zimp Aimp Zeff']);
0074     %
0075     zZi = [1,1,1,Zimp];
0076     zmi = [1,2,3,Aimp];
0077     fi = [frac_H frac_D,0];
0078     %
0079     delete(filename)
0080     %
0081 elseif type == 1.45 ,%.mat, DIII-D profiles
0082     %
0083     eval(['load -mat ',filename,' prho pne pTe pTi frac_H frac_D Zimp Aimp pZeff']);
0084     %
0085     zZi = [1,1,1,Zimp];
0086     zmi = [1,2,3,Aimp];
0087     fi = [frac_H frac_D,0];
0088     %
0089     delete(filename)
0090     %
0091 elseif type == 1.5 ,%.mat, DEMO profile (Artaud)
0092     %
0093     load(filename,'-mat');%prof data save in equil_prof structure
0094     %
0095     pspsin = (equil_prof.psi - equil_prof.psi(1))/(equil_prof.psi(end) - equil_prof.psi(1));
0096     pTe = equil_prof.pTe/1000;%keV
0097     pne = equil_prof.pne;
0098     pzTi = equil_prof.pzTi/1000;%keV
0099     pzni = equil_prof.pzni;
0100     zZi = equil_prof.zZi;
0101     zmi = equil_prof.zmi;
0102     %
0103 elseif type == 1.6 ,%.mat, COMPASS
0104     %
0105     eval(['load -mat ',filename,' pspsin pne pTe pTi Zeff zZi zmi fi']);
0106     %
0107     delete(filename)
0108     %
0109 elseif type == 2,%ascii - 203040
0110     [prho,pTi,pTe,pni,pne] = textread(filename,'%n%n%n%n%n');
0111 elseif type == 3,%ascii formatted (CMOD)
0112     pne = [];
0113     pTe = [];
0114     fid = fopen(filename);
0115     s = 0;
0116     t = 0;
0117     while isempty(s) | s ~= -1,
0118         s = fgetl(fid);
0119         if isempty(findstr(s,'genray')),
0120             continue
0121         else
0122             t = t + 1;
0123             s = fgetl(fid);
0124             s = fgetl(fid);
0125             while isempty(s) | s ~= -1,
0126                 d = sscanf(s,'%f');
0127                 if isempty(d),
0128                     break
0129                 elseif t == 1,
0130                     pne = [pne,d(1)];
0131                 elseif t == 2, 
0132                     pTe = [pTe,d(1)];                    
0133                 end
0134                 s = fgetl(fid);
0135             end
0136         end
0137     end
0138     fclose(fid);
0139     clear fid s t d
0140     %
0141     pne = pne*1e19;
0142     pTi = pTe*2/3;
0143     pspsinT = linspace(0,1,length(pne));
0144     %
0145 elseif type == 3.1,%ascii formatted (CMOD)
0146     pne = [];
0147     pTe = [];
0148     fid = fopen(filename);
0149     s = 0;
0150     while isempty(s) | s ~= -1,
0151         s = fgetl(fid);
0152         if isempty(findstr(s,'Density')),
0153             continue
0154         else
0155             s = fgetl(fid);
0156             s = fgetl(fid);
0157             while ~isempty(s),
0158                 d = sscanf(s,'%f').';
0159                 pne = [pne,d];
0160                 s = fgetl(fid);
0161             end
0162             break
0163         end
0164     end
0165     while isempty(s) | s ~= -1,
0166         s = fgetl(fid);
0167         if isempty(findstr(s,'genray')),
0168             continue
0169         else
0170             s = fgetl(fid);
0171             s = fgetl(fid);
0172             while s ~= -1,
0173                 d = sscanf(s,'%f');
0174                 pTe = [pTe,d(1)];                    
0175                 s = fgetl(fid);
0176             end
0177             break
0178         end
0179     end
0180     fclose(fid);
0181     clear fid s d
0182     %
0183     pTi = pTe*2/3;
0184     pspsinT = linspace(0,1,length(pne));
0185     %
0186 elseif type == 3.2,%ascii formatted (CMOD)
0187     pne = [];
0188     pTe = [];
0189     pTi = [];
0190     fid = fopen(filename);
0191     s = 0;
0192     while isempty(s) | s ~= -1,
0193         s = fgetl(fid);
0194         if isempty(findstr(s,'Density')),
0195             continue
0196         else
0197             s = fgetl(fid);
0198             s = fgetl(fid);
0199             while ~isempty(s),
0200                 d = sscanf(s,'%f').';
0201                 pne = [pne,d];
0202                 s = fgetl(fid);
0203             end
0204             break
0205         end
0206     end
0207     while isempty(s) | s ~= -1,
0208         s = fgetl(fid);
0209         if isempty(findstr(s,'Electron')),
0210             continue
0211         else
0212             s = fgetl(fid);
0213             s = fgetl(fid);
0214             while ~isempty(s),
0215                 d = sscanf(s,'%f').';
0216                 pTe = [pTe,d];
0217                 s = fgetl(fid);
0218             end
0219             break
0220         end
0221     end
0222      while isempty(s) | s ~= -1,
0223         s = fgetl(fid);
0224         if isempty(findstr(s,'Ion')),
0225             continue
0226         else
0227             s = fgetl(fid);
0228             s = fgetl(fid);
0229             while ~isempty(s),
0230                 d = sscanf(s,'%f').';
0231                 pTi = [pTi,d];
0232                 s = fgetl(fid);
0233             end
0234             break
0235         end
0236     end   
0237     fclose(fid);
0238     clear fid s d
0239     %
0240     pne = pne*1e6;%m-3
0241     pspsinT = linspace(0,1,length(pne));
0242     %
0243     zZi = [1,1,1,4];
0244     zmi = [1,2,3,9];
0245     fi = [0,1,0];
0246     %
0247 elseif type == 4,%ascii formatted (ITER)
0248     pspsinT = [];
0249     pne = [];
0250     pTe = [];
0251     pTi = [];
0252     pni = [];
0253     zZi = [];
0254     zmi = [];
0255     fid = fopen(filename);
0256     s = fgetl(fid);
0257     while ~isempty(s),
0258         d = sscanf(s,'%f %f %f %f');
0259         pspsinT = [pspsinT,d(1)];        
0260         pne = [pne,d(2)];
0261         pTe = [pTe,d(3)];
0262         pTi = [pTi,d(4)];
0263         s = fgetl(fid);
0264     end
0265     s = fgetl(fid);
0266     while s ~= -1,
0267         d = sscanf(s,'%f %f %f');
0268         zZi = [zZi,d(1)];
0269         zmi = [zmi,d(2)];
0270         pzni = [pni;pne*d(3)];
0271         s = fgetl(fid);
0272     end
0273     fclose(fid);
0274     clear fid s d
0275     %
0276     pne = pne*1e19;
0277     pzni = pni*1e19;
0278     pTe = pTe;
0279     pTi = pTi;
0280     %
0281     clear pzni % ion densities recalculated from Zeff
0282     fi = [0,0.5,0.5];
0283 elseif type == 5,%ascii formatted (MST)
0284     prho = [];
0285     pspsin = [];
0286     pspsinT = [];
0287     pne = [];
0288     pTe = [];
0289     fid = fopen(filename);
0290     for ii = 1:sscanf(fgetl(fid),'%f'),
0291         d = sscanf(fgetl(fid),'%f %f %f');     
0292         prho = [prho,d(1)];
0293         pne = [pne,d(2)];
0294         pTe = [pTe,d(3)];% in eV units
0295     end
0296     fclose(fid);
0297 %
0298     prho = prho/max(prho);
0299     pTe = pTe/1e3;%keV
0300     zZi = [1,1,1];%No impurity
0301     zmi = [1,2,3];
0302     fi = [0,1,0];%Deuterium plasma
0303     pTi = pTe/2.0;%As done in RANT3D/AORSA
0304 elseif type == 6,%ascii formatted (FTU)
0305     prho = [];
0306     pspsin = [];
0307     pspsinT = [];
0308     pne = [];
0309     pTe = [];
0310     fid = fopen(filename);
0311     s = fgetl(fid);
0312     dummy = sscanf(s,'%s');
0313     s = fgetl(fid);
0314     while s ~= -1,
0315         d = sscanf(s,'%f %f %f');     
0316         pspsin = [pspsin,d(1)];
0317         pne = [pne,d(2)];
0318         pTe = [pTe,d(3)];% in eV units
0319         s = fgetl(fid);        
0320     end
0321     fclose(fid);
0322 %
0323     zZi = [1,1,1];%No impurity
0324     zmi = [1,2,3];
0325     fi = [0,1,0];%Deuterium plasma
0326     pTi = pTe/2.0;%As done in RANT3D/AORSA
0327 elseif type == 7,%ascii formatted (ASCOT)
0328     prho = [];
0329     pspsin = [];
0330     pspsinT = [];
0331     pne = [];
0332     pTe = [];
0333     fid = fopen(filename);
0334     s = fgetl(fid);
0335     dummy = sscanf(s,'%s%c');s = fgetl(fid);
0336     dummy = sscanf(s,'%s%c');s = fgetl(fid);
0337     dummy = sscanf(s,'%s%c');s = fgetl(fid);
0338     d = sscanf(s,'%d %d %s');s = fgetl(fid);nrho = d(1);nion = d(2);
0339     d = sscanf(s,'%d %d %d %d %d %d');s = fgetl(fid);zZi = d';
0340     d = sscanf(s,'%d %d %d %d %d %d');s = fgetl(fid);zmi = d';
0341     dummy = sscanf(s,'%s%c');s = fgetl(fid);
0342     dummy = sscanf(s,'%s%c');s = fgetl(fid);
0343     %
0344     if length(find(zZi==1)) == 1,
0345         pzni = zeros(nion+2,nrho);
0346         pzTi = zeros(nion+2,nrho);
0347     elseif length(find(zZi==1)) == 2,
0348         pzni = zeros(nion+1,nrho);
0349         pzTi = zeros(nion+1,nrho);
0350     else length(find(zZi==1)) == 3,
0351         pzni = zeros(nion,nrho);
0352         pzTi = zeros(nion,nrho);
0353     end
0354     %
0355     j = 1;
0356     %
0357     while s ~= -1,
0358         d = sscanf(s,'%f %f %f %f %f %f %f %f %f %f %f');     
0359         prho = [prho,d(1)];
0360         pTe = [pTe,d(2)/1000];% in eV units
0361         pne = [pne,d(3)];
0362         %
0363         for iion = 1:nion,
0364             if length(find(zZi==1)) == 3,
0365                 pzni(iion,j) = d(5+iion);
0366                 pzTi(iion,j) = d(5)/1000;% in eV units
0367             elseif length(find(zZi==1)) == 2,
0368                 pzni(iion+1,j) = d(5+iion);
0369                 pzTi(iion+1,j) = d(5)/1000;% in eV units
0370             elseif length(find(zZi==1)) == 1,
0371                 pzni(iion+2,j) = d(5+iion);
0372                 pzTi(iion+2,j) = d(5)/1000;% in eV units
0373             end
0374         end
0375         %
0376         j = j+1;
0377         %
0378         s = fgetl(fid);        
0379     end
0380     %
0381     fclose(fid);
0382     %
0383     if length(find(zZi==1)) == 1,
0384         zZi = [1,1,zZi];
0385         zmi = [1,2,zmi];
0386     elseif length(find(zZi==1)) == 2,
0387         zZi = [1,zZi];
0388         zmi = [1,zmi];
0389     end
0390     %
0391     nion = length(zZi);
0392 end
0393 %
0394 if type == -2 || (type < 4 && type >= 0 && type ~= 1.3 && type ~= 1.4 && type ~= 1.45 && type ~= 1.5 && type ~= 1.6 && type ~= 3.2) ,
0395     %
0396     zZi = [1,1,1,6];
0397     zmi = [1,2,3,12];
0398     fi = [0,1,0];
0399     %
0400 end
0401 %
0402 nz = length(zZi);
0403 %
0404 if type == 0 || type == 2,
0405     %
0406     pTe = pTe/1e3;
0407     pTi = pTi/1e3;
0408     %
0409     Zimp = zZi(4);
0410     %
0411     pnimp = (pne-pni)/(Zimp-1);
0412     pnH = pni-pnimp;
0413     %
0414     pzni = [fi'*pnH;pnimp];
0415     pzTi = ones(nz,1)*pTi;
0416     %
0417 elseif type == 0.1,
0418     %
0419     fi = [0,0.5,0.5];
0420     %
0421     psi = psi.';
0422     pTe = pTe.';
0423     pTi = pTi.';
0424     pne = pne.'*1e19;
0425     pZeff = pZeff.';
0426     %
0427     Zimp = zZi(end);
0428     %
0429     pnH = pne.*(Zimp - pZeff)/(Zimp - 1);
0430     pnZ = pne.*(pZeff - 1)/Zimp/(Zimp - 1);
0431     pzni = [fi.'*pnH;pnZ];
0432     pzTi = ones(nz,1)*pTi;
0433     %
0434     pspsin = sqrt(psi);
0435     %
0436 elseif type == 0.2,% note pnH, pnalpha and pnimp not accounted for because impurity species unknown
0437     %
0438     fi = [0,0.5,0.5];
0439 %    zZi = [1,1,1,2,4,6];
0440 %    zmi = [1,2,3,4,9,12];
0441     %
0442     pspsinT = pspsinT.';
0443     pTe = pTe.';
0444     pTi = pTi.';
0445     pne = pne.'*1e19;
0446 %    pnH = pnH.'*1e19;
0447 %    pnalpha = pnalpha.'*1e19;
0448 %    pnimp = pnimp.'*1e19;
0449     clear pnH pnalpha pnimp
0450     pZeff = pZeff.';
0451     %
0452     ZH = zZi(3);
0453     Zalpha = zZi(4);
0454     Zimp1 = zZi(end-1);
0455     Zimp2 = zZi(end);
0456     %
0457 %    pnZimp1 = ((pZeff - Zimp2).*pne - (ZH - Zimp2)*ZH*pnH - (Zalpha - Zimp2)*Zalpha*pnalpha)/Zimp1/(Zimp1 - Zimp2);
0458 %    pnZimp2 = ((pZeff - Zimp1).*pne - (ZH - Zimp1)*ZH*pnH - (Zalpha - Zimp1)*Zalpha*pnalpha)/Zimp2/(Zimp2 - Zimp1);
0459 %    pzni = [fi.'*pnH;pnalpha;pnZimp1;pnZimp2];
0460     Zimp = zZi(end);
0461     pnH = pne.*(Zimp - pZeff)/(Zimp - 1);
0462     pnZ = pne.*(pZeff - 1)/Zimp/(Zimp - 1);
0463     pzni = [fi.'*pnH;pnZ];
0464     pzTi = ones(nz,1)*pTi;
0465     %
0466 elseif type == -2 || type == 1,%TCV
0467     %
0468     pTe = pTe/1e3;
0469     pTi = pTi/1e3;
0470     %
0471     pspsin = prho;% rho definition in TCV
0472     prho = [];
0473     %
0474     pnH = pne.*(zZi(4) - pZeff)/(zZi(4) - 1);
0475     pnZ = pne.*(pZeff - 1)/zZi(4)/(zZi(4) - 1);
0476     pzni = [fi.'*pnH;pnZ];
0477     pzTi = ones(nz,1)*pTi;
0478 end
0479 if type == 1.1 || type == 1.2 || type == 1.4 || type == 1.6 || type>= 3 && type ~=5 && type ~=6 && type ~=7,
0480     %
0481     pnH = pne*(zZi(4) - Zeff)/(zZi(4) - 1);
0482     pnZ = pne*(Zeff - 1)/zZi(4)/(zZi(4) - 1);
0483     pzni = [fi.'*pnH;pnZ];
0484     pzTi = ones(nz,1)*pTi;
0485 elseif type == 1.45,
0486     pnH = pne.*(zZi(4) - pZeff)/(zZi(4) - 1);
0487     pnZ = pne.*(pZeff - 1)/zZi(4)/(zZi(4) - 1);
0488     pzni = [fi.'*pnH;pnZ];
0489     pzTi = ones(nz,1)*pTi;
0490 elseif type == 0 || type == 2,
0491     pnH = pne;
0492     pzni = [fi.'*pnH];
0493     pzTi = ones(nz,1)*pTi;
0494 end
0495 %
0496 if display_mode >= 2,
0497     %
0498     style = '-';
0499     marker = '+';
0500     color = 'r';
0501     width = 2;
0502     siz = 20;
0503     red = 0.9;
0504     lspace = 0.7;
0505     bspace = 0.6;
0506     tit = strrep(dispname,'_','-');
0507     %
0508     if ~isempty(prho),
0509         xlab = '\rho';
0510         pdisp = prho;
0511     elseif ~isempty(pspsin)
0512         xlab = '\psi_n^{1/2}';
0513         pdisp = pspsin;
0514     elseif ~isempty(pspsinT)
0515         xlab = '\psi_{nT}^{1/2}';
0516         pdisp = pspsinT;
0517     end        
0518     %
0519     figure,clf
0520     %
0521     ylab = 'T_e (keV)';
0522     xlim = [0,1];
0523     ylim = NaN;
0524     %
0525     graph1D_jd(pdisp,pTe,0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker,color,width,siz,gca,red,lspace,bspace);
0526     set(gca,'XTick',[0:0.2:1]);
0527     set(gca,'ylim',[0,max(get(gca,'ylim'))])
0528     print_jd(p_opt,[savename,'_Te_rho']);
0529     %
0530     figure,clf
0531     %
0532     ylab = 'n_e (\times 10^{19} m^{-3})';
0533     %
0534     graph1D_jd(pdisp,pne/1e19,0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker,color,width,siz,gca,red,lspace,bspace);
0535     set(gca,'XTick',[0:0.2:1]);
0536     set(gca,'ylim',[0,max(get(gca,'ylim'))])
0537     print_jd(p_opt,[savename,'_ne_rho']);
0538     %
0539     figure,clf
0540     %
0541     ylab = 'T_i (keV)';
0542     %
0543     graph1D_jd(pdisp,pzTi,0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker,color,width,siz,gca,red,lspace,bspace);
0544     set(gca,'XTick',[0:0.2:1]);
0545     set(gca,'ylim',[0,max(get(gca,'ylim'))])
0546     print_jd(p_opt,[savename,'_Ti_rho']);
0547     %
0548     figure,clf
0549     %
0550     ylab = 'n_i (\times 10^{19} m^{-3})';
0551     %
0552     graph1D_jd(pdisp,pzni/1e19,0,0,xlab,ylab,tit,NaN,xlim,ylim,style,marker,color,width,siz,gca,red,lspace,bspace);
0553     set(gca,'XTick',[0:0.2:1]);
0554     set(gca,'ylim',[0,max(get(gca,'ylim'))])
0555     print_jd(p_opt,[savename,'_ni_rho']);
0556     %
0557  end
0558 
0559 
0560 
0561 
0562 
0563 
0564 
0565