load_externaldata_JET

 This function loads data from JET database using zjet

0001 function external = load_externaldata_JET(workdir,shotnum,shotime,t1,t2,opt_gui)
0002 %
0003 % This function loads data from JET database using zjet
0004 %
0005 external = '';
0006 % basestr = 'JET';
0007 %
0008 % -------------------------------------------------------------------------
0009 %
0010 % addpath('/luke/Projet_Cronos/acces/jet/interface')
0011 % addpath('/luke/Projet_Cronos/')
0012 % zineb_path
0013 % %
0014 % % script de demarrage de l'interface de preparation des donnees JET
0015 % %
0016 % setappdata(0,'langue_cronos','anglais');
0017 % %
0018 % if ~exist('zuijetdata','file'),
0019 %     disp('-----> Access to the JET database is not granted.')
0020 %     disp('-----> You need the zuijetdata.m file in the matlab path.');
0021 %     return;
0022 % end
0023 % %
0024 % zuijetdata(shotnum);
0025 %
0026 % -------------------------------------------------------------------------
0027 %
0028 % filestr = [shotnum,'_',shotime,'_',num2str(t1),'_',num2str(t2)];
0029 % datafile = ['data_for_LUKE_',filestr,'.mat'];
0030 % %
0031 % [datafile,datapath] = igetfile_jd(opt_gui,datafile,'Please select the JET data file',workdir);
0032 % %
0033 % if datafile == 0,
0034 %     return
0035 % end
0036 % %
0037 % load([datapath,datafile],'x','y','Rp','Zp','xypsin','psia','xBt',...
0038 %     'xne','xTe','Zeff');
0039 %
0040 % -------------------------------------------------------------------------
0041 %
0042 addpath('/luke/Projet_Cronos/acces/jet/bilan')
0043 addpath('/luke/Projet_Cronos/acces/jet/bilan/zwingmann')
0044 %
0045 shotimenum = str2num(shotime);
0046 %
0047 datadir = igetdir_jd(opt_gui,'Please select the JET data directory',workdir);
0048 %
0049 efitstr = [datadir,'/efit',shotnum,'.mat'];
0050 tempstr = [datadir,'/temp',shotnum,'.mat'];
0051 profstr = [datadir,'/prof',shotnum,'.mat'];
0052 nexstr = [datadir,'/nex',shotnum,'.mat'];
0053 texstr = [datadir,'/tex',shotnum,'.mat'];
0054 %
0055 if exist(efitstr,'file') && exist(tempstr,'file') && exist(profstr,'file'),
0056     %
0057     load(tempstr,'emap','fbnd','faxs','tzeff','zefv','zefh');
0058     load(efitstr,'tefit','psix','qx');
0059     load(profstr,'rbnd','zbnd','rth','tth','neth','teth');
0060     %
0061     if exist(nexstr,'file')
0062         load(nexstr,'nex');
0063     else
0064         nex = [];
0065     end
0066     %
0067     if exist(texstr,'file')
0068         load(texstr,'tex');
0069     else
0070         tex = [];
0071     end
0072     %
0073 else
0074     error(['equilibrium data not found for shot ',shotnum]);
0075 end
0076 %
0077 % Magnetic Equilibrium
0078 %
0079 timediff = tefit - shotimenum;
0080 kf = find(abs(timediff) == min(abs(timediff)),1,'first');
0081 psia = fbnd(kf);
0082 psip_old = faxs(kf);
0083 psix = psix(kf,:);
0084 pq = qx(kf,:);
0085 %
0086 np = 100;
0087 nz = 100;
0088 %
0089 R2 = linspace(min(rbnd(kf,:))*0.9,max(rbnd(kf,:))*1.1,np);
0090 Z2 = linspace(min(zbnd(kf,:))*1.1,max(zbnd(kf,:))*1.1,nz);
0091 [R,Z] = meshgrid(R2,Z2);
0092 %
0093 [psi,Br,Bz,Bt,err] = jetequi3(R,Z,emap,tefit(kf));%+40
0094 %
0095 xypsi = psi.';
0096 %xyBr = Br.';
0097 %xyBz = Bz.';
0098 xBt = -Bt;% negative sign for phi orientation
0099 %
0100 R = R(1,:);
0101 Z = Z(:,1).';
0102 %
0103 mask = find(all(xypsi.' == repmat(xypsi(:,1).',[nz,1])));
0104 %
0105 R(mask) = [];
0106 xypsi(mask,:)=  [];
0107 %xyBr(mask,:)=  [];
0108 %xyBz(mask,:)=  [];
0109 xBt(mask)=  [];
0110 %
0111 sigma = sign(psia);
0112 Ns = 101;
0113 %
0114 psip = min(min(sigma*xypsi));
0115 [iR,iZ] = find(sigma*xypsi == psip);
0116 %
0117 if length(iR) ~= 1,
0118     error('unable to locate minimum in xypsi data');
0119 end
0120 %
0121 ss = linspace(0,1,Ns);
0122 sR = R(iR-1) + (R(iR+1)-R(iR-1))*ss;
0123 sZ = Z(iZ-1) + (Z(iZ+1)-Z(iZ-1))*ss;
0124 %
0125 sspsi = interp2(R,Z,xypsi.',sR.',sZ,'spline').';
0126 psip = sigma*min(min(sigma*sspsi));
0127 [isR,isZ] = find(sspsi == psip);
0128 %
0129 if length(isR) ~= 1,
0130     error('unable to locate minimum in xypsi data');
0131 end
0132 %
0133 Rp = sR(isR);
0134 Zp = sZ(isZ);
0135 %
0136 disp(['Old value of psip: ',num2str(psip_old),', new value: ',num2str(psip)])
0137 disp(['Rp: ',num2str(Rp)])
0138 disp(['Zp: ',num2str(Zp)])
0139 %
0140 x = R - Rp;
0141 y = Z - Zp;
0142 xypsi = xypsi - psip;
0143 psia = psia - psip;
0144 psix = psix - psip;
0145 %
0146 xypsin = xypsi/psia;
0147 psin = psix/psia;
0148 psin(1) = 0;
0149 %
0150 % Profiles
0151 %
0152 xth = rth - Rp;
0153 method = 'linear';
0154 %
0155 nt = length(tth);
0156 %
0157 tth = tth.';
0158 %
0159 tth_S = [tth(1),(tth(1:end-1) + tth(2:end))/2,tth(end)];
0160 dtth = min([tth_S(2:end);t2*ones(1,nt)]) - max([tth_S(1:end-1);t1*ones(1,nt)]);
0161 dtth(dtth < 0) = 0;
0162 rdtth = repmat(dtth.',[1,length(rth)]);
0163 pdtth = repmat(dtth.',[1,length(psin)]);
0164 %
0165 neth = sum(neth.*rdtth)/(t2-t1);
0166 teth = sum(teth.*rdtth)/(t2-t1);
0167 if ~isempty(nex),
0168     pne = sum(nex.*pdtth)/(t2-t1);
0169 else
0170     pne = [];
0171 end
0172 if ~isempty(tex),
0173     pTe = sum(tex.*pdtth)/(t2-t1);
0174 else
0175     pTe = [];
0176 end
0177 %
0178 xne = interp1(xth,neth,x,method);
0179 xTe = interp1(xth,teth,x,method);
0180 Zeff = interp1(tzeff,(zefh + zefv)/2,shotimenum);
0181 %
0182 filestr = [shotnum,'_',shotime,'_',num2str(t1),'_',num2str(t2),'.mat'];
0183 %
0184 % Data saved for irunluke
0185 %
0186 external.equil.shotnum = shotnum;
0187 external.equil.shotime = shotime;
0188 %
0189 external.equil.magnetic.x = x;
0190 external.equil.magnetic.y = y;
0191 external.equil.magnetic.Rp = Rp;
0192 external.equil.magnetic.Zp = Zp;
0193 external.equil.magnetic.xypsin = xypsin;
0194 external.equil.magnetic.psia = psia;
0195 external.equil.magnetic.xBt = xBt;
0196 external.equil.magnetic.psin = psin;
0197 external.equil.magnetic.pq = pq;
0198 %
0199 external.equil.prof.xne = xne;
0200 external.equil.prof.xTe = xTe;
0201 external.equil.prof.Zeff = Zeff;
0202 external.equil.prof.psin = psin;
0203 external.equil.prof.pne = pne;
0204 external.equil.prof.pTe = pTe;
0205 %
0206 % external.ohm. = ;
0207 % external.wave. = ;
0208 %
0209 external.id = ['JET_',filestr];
0210