/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 21.05.2 ] */
/* [wxMaxima: hide output   ] *//* [wxMaxima: input   start ] */
kill(all);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
/*　ユーザ定義関数：対角行列化　*/
Diag(m):=ident(length(m))*m$;

/*　ユーザ定義関数：クリストッフェルの3指標記号（1種）*/
C3f(gxx,X):=makelist(makelist(makelist(
             (1/2)*(diff(gxx[i,j],X[k])+diff(gxx[i,k],X[j])-diff(gxx[j,k],X[i])),
           k,1,length(gxx)),j,1,length(gxx)),i,1,length(gxx))$

/*　ユーザ定義関数：クリストッフェルの3指標記号（2種）*/
C3s(gxx,X):=makelist(makelist(makelist(
             sum(invert(gxx)[i,l]*C3f(gxx,X)[l][j][k],l,1,length(gxx)),
           k,1,length(gxx)),j,1,length(gxx)),i,1,length(gxx))$

/*　ユーザ定義関数：リーマンク・リストッフェルの曲率テンソル(計量指定形）　*/
R4_g(gxx,X):=makelist(makelist(makelist(makelist(
             diff(C3s(gxx,X)[h][k][i],X[j])-
             diff(C3s(gxx,X)[h][j][i],X[k])+
             sum(C3s(gxx,X)[h][j][l]*C3s(gxx,X)[l][k][i],l,1,length(gxx))-
             sum(C3s(gxx,X)[h][k][l]*C3s(gxx,X)[l][j][i],l,1,length(gxx)),
           k,1,length(gxx)),j,1,length(gxx)),i,1,length(gxx)),h,1,length(gxx))$
    
/*　ユーザ定義関数：リーマンク・リストッフェルの曲率テンソル(接続指定形）　*/
R4(C3s,X):=makelist(makelist(makelist(makelist(
             diff(C3s[h][k][i],X[j])-
             diff(C3s[h][j][i],X[k])+
             sum(C3s[h][j][l]*C3s[l][k][i],l,1,length(C3s))-
             sum(C3s[h][k][l]*C3s[l][j][i],l,1,length(C3s)),
           k,1,length(C3s)),j,1,length(C3s)),i,1,length(C3s)),h,1,length(C3s))$
        
/*　ユーザ定義関数：リッチテンソル(曲率指定形）　*/
R2(R4):=makelist(makelist(
              sum(R4[h][i][j][h],h,1,length(R4)),
            j,1,length(R4)),i,1,length(R4))$
            
/*　ユーザ定義関数：スカラー曲率(曲率指定形）　*/
R1(R2,gxx):=sum(sum(invert(gxx)[i][j]*R2[i][j],i,1,length(R2)),j,1,length(R2))$
                
/*　ユーザ定義関数：アインシュタインテンソル　*/
Gxx(R4,gxx):=R2(R4)-(1/2)*gxx*R1(R2(R4),gxx)$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
/*　ユーザ定義関数：自明解フィルター　*/
  TrivFilter2d( In):=block([sz,Out],
    sz:length(In), Out:[],
    for i:1 thru sz do(for j:1 thru sz do(
      if i>=j then if is(trigsimp(trigreduce(In[i][j]))#-trigsimp(trigreduce(In[j][i]) ))  then
      /* Not equal case or Unknown case... */
            Out:endcons((In[i][j]+In[j][i])=0, Out)
    )),
    return(Out)
 )$

/*　ユーザ定義関数：測地線方程式　*/
GeoCvEq(gxx,X,tau):=block([Xtau,v4D,Res],
  Xtau:makelist(X[i](tau),i,length(X)),
  v4D:diff(Xtau,tau),
  for i:1 thru length(X) do(gxx:subst(Xtau[i],X[i],gxx)),
  Res:diff(gxx.v4D,tau)-(1/2)*makelist(v4D.diff(gxx,Xtau[i]).v4D,i,length(X)),
  return(makelist(Res[i][1]=0,i,length(Res)))
)$


/*　ユーザ定義関数：キリングベクトル確認　*/
KillingCheck(v,gxx,C3s,X):=block([sz,Kv],
  sz:length(v),   
  Kv:makelist(makelist(diff((gxx.v)[i,1],X[j]),j,sz),i,sz)
      -sum(C3s[s]*(gxx.v)[s,1],s,1,sz),
  return(TrivFilter2d(Kv))
)$

/*    RuleList[i][2] ->  RuleList[i][1] */
ListSubst(RuleList, Expr):=block([work],
  work:Expr,
  for i:1 thru length(RuleList) do(work: subst(RuleList[i][1],RuleList[i][2],work)),
  return(ev(work,nouns))
)$

/*    RuleList[i][2] ->  RuleList[i][1] */
ListRatsubst(RuleList, Expr):=block([work],
  work:Expr,
  for i:1 thru length(RuleList) do(work: ratsubst(RuleList[i][1],RuleList[i][2],work)),
  return(ev(work,nouns))
)$

Dimensions(List):=block([temp],
  temp:[length(List)],
  for ix:1 while is(length(List)#length(flatten(List))) do(
    temp:endcons(length(List[1]),temp),
    List:List[1]
  ),
  return(temp)
)$

PrintNonZero(List,ID):=block([ix],
  for ix:1 thru length(List) do(
    if is(length(List)#length(flatten(List))) then
      PrintNonZero(List[ix],concat(ID,ix-1))
    else
      if is(List[ix]#0)then
        print(concat(ID,(ix-1)),"=",List[ix])
  )
)$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
s;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
swzRuleList:[ g00=lambda([p1] ,(1-2*m/p1)), 
              g11=lambda([p1] ,(-1/(1-2*m/p1)))];
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
kt:[1,0,0,0];
kz:[0,0,0,-1];
kx:[0,0,-sin(phi),-cos(phi)*cot(theta)];
ky:[0,0,cos(phi),-sin(phi)*cot(theta)];
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
/*map(lambda([p1],KillingCheck(p1,swzgxx,swzC3s,X)),[kt,kz,kx,ky]);*/;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
V:'[diff(t,s),diff(r,s),diff(theta,s),diff(phi,s)];
kConst:[en=V.swzgxx.kt, jx=V.swzgxx.kx, jy=V.swzgxx.ky, jz=V.swzgxx.kz];
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
swzX:[x0,r,theta,phi]$
swzgxx:Diag(
    [g00(r),-g11(r),-r^2,-r^2*sin(theta)^2]);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
swzC3s:C3s(swzgxx,swzX),trigreduce$
PrintNonZero(swzC3s,"C3s_")$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
swzR4:R4(swzC3s,swzX),ratsimp,ratreduce$
swzR4:swzR4,expand$
PrintNonZero(swzR4,"R4_")$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
swzR2:R2(swzR4)$
PrintNonZero(swzR2,"R2_");
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
swzGxx:(trigsimp(expand(Gxx(swzR4,swzgxx))))$
swzGxxUL:ratsimp((expand(args(invert(swzgxx).swzGxx))))$
halfangles:false$ trigexpand:true$
swzGxxUL:trigsimp(swzGxxUL),factor$
PrintNonZero(swzGxxUL,"Gul_");
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
g11eqL:ode2(swzGxxUL[1][1]=0,g11(r),r);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
g11eq:solve(logcontract(g11eqL),g11(r));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
swzg11F:subst(a,%e^%c,g11eq)[1];
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
g00eq:ode2(swzGxxUL[1][1]=swzGxxUL[2][2],g00(r),r);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
swzg00F:subst([%c=b,swzg11F],g00eq);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
swzgxxF:subst([swzg00F,swzg11F],swzgxx);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
ntXc:[x0,x1,x2,x3];
ntXf:[X0,X1,X2,X3];
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
ruleXcXf:[X0=x0,X1=x1-(h-(1/2)*g*x0^2/c^2),X2=x2,
        X3=x3];
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
ntXcf:subst(ruleXcXf,ntXf);
Jh:apply(matrix,makelist(makelist(diff(ntXcf[ix],ntXc[iy]),iy,4),ix,4))$; 
eta:ident(4)*[1,-1,-1,-1]$;
ntgxx:transpose(Jh).eta.Jh;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
ntPhi:g*(x1-h);
Phieq:Phi=subst(x1=h-(1/2)*g*x0^2/c^2 ,ntPhi);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
x0x0eq:solve(Phieq,x0^2)[1];
ntg00F:subst([x0x0eq,Phi=-G*M/r],ntgxx[1][1] );
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
ntg00F=expand(swzgxxF[1][1]);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
aeq:subst(b=1,solve(ntg00F=swzgxxF[1][1],a));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
swzXFr:subst(x0=c*t,swzX);
swzgxxFr:subst([b=1,aeq[1]],swzgxxF),expand;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
const:[M=2.0E33,G=6.7E-14,c=3.0E8]$
redshift:subst(const,solve(swzgxxFr[1][1]=0,r));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
horizon:subst(const,solve(1/swzgxxFr[2][2]=0,r));
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
swzgxxFr[2];
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
wxplot3d( [100,subst(horizon,r),subst(redshift,r),
  [theta, 0, %pi], [phi, 0, 4*%pi/4]],
  [transform_xy, spherical_to_xyz],  
  [grid,20,20],[legend,false],[same_xyz,true])$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
kRuleList:[
  g00=lambda([p1,p2], 
    (p1^2-2*m*p1+a^2)/( p1^2+a^2*cos(p2)^2)
    -a^2*sin(p2)^2/( p1^2+a^2*cos(p2)^2)),
  g03=lambda([p1,p2],
    -(p1^2-2*m*p1+a^2)/( p1^2+a^2*cos(p2)^2)
    *a*sin(p2)^2
    +sin(p2)^2/( p1^2+a^2*cos(p2)^2)*(p1^2+a^2)*a),
  g11=lambda([p1,p2],
    -(p1^2+a^2*cos(p2)^2)/(p1^2-2*m*p1+a^2) ),
  g22=lambda([p1,p2],
    -(p1^2+a^2*cos(p2)^2) ),
  g33=lambda([p1,p2], 
    (p1^2-2*m*p1+a^2)/( p1^2+a^2*cos(p2)^2)
    *a^2*sin(p2)^4
    -sin(p2)^2/( p1^2+a^2*cos(p2)^2)*(p1^2+a^2)^2)
]$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
X:[t,r,theta,phi];
kgxx:Diag([g00(r,theta),g11(r,theta),g22(r,theta),g33(r,theta)])$
kgxx[4][1]:g03(r,theta)$ kgxx[1][4]:g03(r,theta)$
kgxx;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
kC3s:C3s(kgxx,X)$
kC00:kC3s[1][1],epand,ratsimp;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
ans:subst(kRuleList,kC00)$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
halfangles:false;trigexpand:true;
ans1:trigsimp(trigrat(ratsimp(ans[2]))),factor$
ans2:trigsimp(ans[2]),factor$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
eq1:ratsimp(subst(kRuleList,g00(r,theta)))=0;
redshift:solve(eq1,r);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
eq2:invert(kgxx)[2][2]=0,ratsimp;
horizon:solve(subst(kRuleList,eq2),r);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
fig[3];
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
fH;
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
fig:map( lambda( [p1], subst(p1,r)), 
  append(horizon,redshift) )$
fH: append(subst([a=0.99,m=1],[fig[1],fig[2]]), 
                   [[theta, 0, %pi], [phi, 0, 4*%pi/4]] );
wxplot3d( 
  fH ,[palette, [gradient, red,  yellow]],
  [transform_xy, spherical_to_xyz],  
  [grid,20,20],[legend,false],[same_xyz,true])$
fR: append(subst([a=0.99,m=1],[fig[3],fig[4]]), 
                   [[theta, 0, %pi], [phi, 0, 4*%pi/4]] );
wxplot3d( 
  fR ,[palette, [gradient,  green, blue]],
  [transform_xy, spherical_to_xyz],  
  [grid,20,20],[legend,false],[same_xyz,true])$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
example(wxplot3d);
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
apropos("wxplot3d");
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
wxplot3d( 
  append(
    subst([a=0*m,m=M*G/c^2,M=2.0E33,
           G=6.7E-14,c=3.0E8],fig),
    [[theta, 0, %pi], [phi, 0, 4*%pi/4]]
  ),
  [transform_xy, spherical_to_xyz],  
  [grid,20,20],[legend,false],[same_xyz,true]
)$
/* [wxMaxima: input   end   ] */


/* [wxMaxima: input   start ] */
wxplot3d( 
  append(
    subst([a=0.99*m,m=M*G/c^2,M=2.0E33,
           G=6.7E-14,c=3.0E8],fig),
    [[theta, 0, %pi], [phi, 0, 4*%pi/4]]
  ),
 [transform_xy, spherical_to_xyz],  
 [grid,20,20],[legend,false],[same_xyz,true])$
/* [wxMaxima: input   end   ] */



/* Old versions of Maxima abort on loading files that end in a comment. */
"Created with wxMaxima 21.05.2"$