【高强度聚焦超声模拟器】模拟分层介质中的高强度聚焦超声波束和加热效应（Matlab代码）

function [yplot]=prueba2(xpun,C,R,w)xpun=[ones(1,size(xpun,2));xpun];yplot=zeros(1,size(xpun,2));for i=1:size(xpun,2)    yplot(i)=0;    for n=1:size(C,1)        z=0;        for a=1:size(xpun,1)            z=z+(xpun(a,i)-C(n,a))^2/R^2;        end        yplot(i)=yplot(i)+w(n)*exp(-z);    endendeval(['load ','''',[pathname,Xcalname,'''']])factores=eval(strtok(Xcalname,'.'));eval(['load ','''',[pathname,ycalname,'''']])y_s=eval(strtok(ycalname,'.'));Nresp=size(y_s,2);Nf=size(factores,2);w_s=w_s/sum(w_s);eval(['load ','''',[pathname,netname,'''']])for i=1:Nf    x=factores(:,i);    delta=(max(x)-min(x));    MM=min(x)+max(x);    x=(2*x-MM)/delta;    eval(['x',int2str(i),'=x;'])    eval(['delta',int2str(i),'f=delta;'])    eval(['MM',int2str(i),'f=MM;'])endNresp=size(y_s,2);for i=1:Nresp    y=y_s(:,i);    delta=(max(y)-min(y));    MM=min(y)+max(y);    y=(2*y-MM)/delta;    eval(['y',int2str(i),'=y;'])    eval(['deltay',int2str(i),'=delta;'])    eval(['MMy',int2str(i),'=MM;'])endfor i_r=1:Nresp    Yplot{i_r}=prueba2(xpun',C{i_r},R(i_r),w{i_r});    eval(['Des',int2str(i_r),'=Yplot{',int2str(i_r),'};'])endfor i_r=1:Nresp    if isnan(lows(i_r))==1        eval(['lows(',int2str(i_r),')=(min(Yplot{',int2str(i_r),'}*deltay',int2str(i_r),'+MMy',int2str(i_r),')/2);'])    end    if isnan(highs(i_r))==1        eval(['highs(',int2str(i_r),')=(max(Yplot{',int2str(i_r),'}*deltay',int2str(i_r),'+MMy',int2str(i_r),')/2);'])    end    switch decides(i_r)        case 1            eval(['Ycomp',int2str(i_r),'=(Yplot{',int2str(i_r),'}*deltay',int2str(i_r),'+MMy',int2str(i_r),')/2;'])            eval(['Ytemp=Ycomp',int2str(i_r),';'])            Dtemp=zeros(size(Ytemp));            Dtemp(Ytemp<lows(i_r))=0;            Dtemp(Ytemp>highs(i_r))=0;            Dtemp(Ytemp>=lows(i_r) & Ytemp<=highs(i_r))=...                (Ytemp(Ytemp>=lows(i_r) & Ytemp<=highs(i_r))-lows(i_r))/(highs(i_r)-lows(i_r));            eval(['Des',int2str(i_r),'=Dtemp;'])        case 2            eval(['Ycomp',int2str(i_r),'=(Yplot{',int2str(i_r),'}*deltay',int2str(i_r),'+MMy',int2str(i_r),')/2;'])            eval(['Ytemp=Ycomp',int2str(i_r),';'])            Dtemp=zeros(size(Ytemp));            Dtemp(Ytemp<lows(i_r))=0;            Dtemp(Ytemp>highs(i_r))=0;            Dtemp(Ytemp>=lows(i_r) & Ytemp<=highs(i_r))=...                (highs(i_r)-Ytemp(Ytemp>=lows(i_r) & Ytemp<=highs(i_r)))/(highs(i_r)-lows(i_r));            eval(['Des',int2str(i_r),'=Dtemp;'])        case 3            eval(['Ycomp',int2str(i_r),'=(Yplot{',int2str(i_r),'}*deltay',int2str(i_r),'+MMy',int2str(i_r),')/2;'])            eval(['Ytemp=Ycomp',int2str(i_r),';'])            Dtemp=zeros(size(Ytemp));            Dtemp(Ytemp<lows(i_r))=0;            Dtemp(Ytemp>highs(i_r))=0;            Dtemp(Ytemp>=lows(i_r) & Ytemp<=targets(i_r))=...                (Ytemp(Ytemp>=lows(i_r) & Ytemp<=targets(i_r))-lows(i_r))/(targets(i_r)-lows(i_r));            Dtemp(Ytemp>targets(i_r) & Ytemp<highs(i_r))=...                (highs(i_r)-Ytemp(Ytemp>targets(i_r) & Ytemp<highs(i_r)))/(highs(i_r)-targets(i_r));            eval(['Des',int2str(i_r),'=Dtemp;'])        case 4            eval(['Ycomp',int2str(i_r),'=(Yplot{',int2str(i_r),'}*deltay',int2str(i_r),'+MMy',int2str(i_r),')/2;'])            eval(['Ytemp=Ycomp',int2str(i_r),';'])            Dtemp=ones(size(Ytemp));            Dtemp(Ytemp<lows(i_r))=0;            Dtemp(Ytemp>highs(i_r))=0;            eval(['Des',int2str(i_r),'=Dtemp;'])    end    % Cleaning bad values    eval(['Des',int2str(i_r),'(Des',int2str(i_r),'<0)=0;'])    eval(['Des',int2str(i_r),'(Des',int2str(i_r),'>1)=1;'])endDes=1;for i_r=1:Nresp    eval(['Des=[Des.*(Des',int2str(i_r),'.^w_s(',int2str(i_r),'))];'])end% Restrictions on factors% eval(['load Dmaxmin',int2str(Nf)])xtempmax=linspace(0,1,15)';xtempmin=linspace(1,0,15)';switch Nf    case 2        dmax2=[kron(xtempmax,ones(15,1)),kron(ones(15,1),xtempmax)];        dmin2=[kron(xtempmin,ones(15,1)),kron(ones(15,1),xtempmin)];    case 3        dmax3=[kron(xtempmax,kron(ones(15,1),ones(15,1))),...            kron(ones(15,1),kron(xtempmax,ones(15,1))),...            kron(ones(15,1),kron(ones(15,1),xtempmax))];        dmin3=[kron(xtempmin,kron(ones(15,1),ones(15,1))),...            kron(ones(15,1),kron(xtempmin,ones(15,1))),...            kron(ones(15,1),kron(ones(15,1),xtempmin))];    case 4        dmax4=[kron(kron(xtempmax,kron(ones(15,1),ones(15,1))),ones(15,1)),...            kron(kron(ones(15,1),kron(xtempmax,ones(15,1))),ones(15,1)),...            kron(kron(ones(15,1),kron(ones(15,1),xtempmax)),ones(15,1)),...            kron(kron(ones(15,1),kron(ones(15,1),ones(15,1))),xtempmax)];        dmin4=[kron(kron(xtempmin,kron(ones(15,1),ones(15,1))),ones(15,1)),...            kron(kron(ones(15,1),kron(xtempmin,ones(15,1))),ones(15,1)),...            kron(kron(ones(15,1),kron(ones(15,1),xtempmin)),ones(15,1)),...            kron(kron(ones(15,1),kron(ones(15,1),ones(15,1))),xtempmin)];endfor i_f=1:Nf    if isnan(lowsf(i_f))==1        eval(['lowsf(',int2str(i_f),')=min(factores(:,',int2str(i_f),'));'])    end    if isnan(highsf(i_f))==1        eval(['highsf(',int2str(i_f),')=max(factores(:,',int2str(i_f),'));'])    end    switch decidesf(i_f)        case 1            eval(['Des=Des.*dmax',int2str(Nf),'(:,i_f)'';'])        case 2            eval(['Des=Des.*dmin',int2str(Nf),'(:,i_f)'';'])        case 3            dtargf=zeros(size(Des));            eval(['xpunf=xpun(:,',int2str(i_f),');'])            eval(['xpunf=(xpunf*delta',int2str(i_f),'f+MM',int2str(i_f),'f)/2;'])            dtargf(xpunf>lowsf(i_f) & xpunf<=targetsf(i_f))=(xpunf(xpunf>lowsf(i_f) & xpunf<=targetsf(i_f))-lowsf(i_f))/(targetsf(i_f)-lowsf(i_f));            dtargf(xpunf>targetsf(i_f) & xpunf<highsf(i_f))=(highsf(i_f)-xpunf(xpunf>targetsf(i_f) & xpunf<highsf(i_f)))/(highsf(i_f)-targetsf(i_f));            Des=Des.*dtargf;        case 4            if lowsf(i_f)>min(factores(:,i_f)) && highsf(i_f)<max(factores(:,i_f))                dtargf=ones(size(Des));                eval(['xpunf=xpun(:,',int2str(i_f),');'])                eval(['xpunf=(xpunf*delta',int2str(i_f),'f+MM',int2str(i_f),'f)/2;'])                dtargf(xpunf<lowsf(i_f))=0;                dtargf(xpunf>highsf(i_f))=0;                Des=Des.*dtargf;            end    endend[MaxDes,maxind]=max(Des);